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A SMOOTHING ALGORITHM USING CUBIC SPLINE FUNCTIONS 


By Robert E. Smith, Jr., Joseph M. Price, and Lona M. Howser 
Langley Research Center 

SUMMARY 

A technique for smoothing arbitrary sets of two-dimensional data is described. 

Both explicit and parametric cubic spline functions are used in a least -squares algorithm 
to obtain a least -squares cubic polynomial spline fit for smoothing data. The primary 
advantage of the least-squares spline technique is that a set of data can be represented by 
a continuous function with continuous first and second derivatives. The technique has 
been programed in FORTRAN IV for the Control Data series 6000 computer systems. 
Interactive computer graphics using a CDC 250 series CRT console and the Control 
Data series 6000 computer graphics system at the Langley Research Center is incor- 
porated into the program to allow a user to interact with the program to aid in obtaining 
a satisfactory solution. A description of the computer program and the procedures for 
using it are included. Examples of applications are presented. 

INTRODUCTION 

Least-squares polynomial curve fitting is a standard smoothing technique for func- 
tionally representing a data set S = The process consists of finding the 

polynomial coefficients which minimize the sum of the squares of the differences between 
the polynomial curve and the values of the dependent variable in the data set. In order 
to obtain a better representation of a data set, polynomials can be fitted to a sequence of 
subsets of the data. The disadvantage of this type of polynomial fit — even though it may 
provide a closer representation of the data - is that there is no guarantee of continuity in 
either the polynomial representation or its derivatives over the entire data set. In con- 
trast to this, applying the least-squares criteria to low-order spline polynomials guar- 
antees continuity in the functional representation and certain of its derivatives while pro- 
ducing a good representation of the data set. 

This concept is discussed by de Boor in references 1 and 2. He notes that the suc- 
cess of using spline functions in smoothing data lies in the proper choice of the "joints." 

In reference 2 he presents an algorithm to choose the number and position of the joints. 

His algorithm is limited as pointed out in reference 2 and is subject to failure on certain 
sets of data. The algorithm and computer program as described in this report emphasizes 



the user interaction with the algorithm to obtain a satisfactory solution. It should be 
noted, however, certain aspects of de Boor's algorithm (ref. 2) would work well in an 
interactive environment to aid the user selection of joints. The basic difference between 
de Boor's fixed joint algorithm (ref. 1) and Smith's algorithm (ref. 3) is the spline func- 
tion representation and use of Lagrange multipliers. All three algorithms would produce 
similar results for data where globally large gradients are not encountered. In this 
report, spline smoothing is extended to data sets which may wander arbitrarily through 
a two-dimensional coordinate space and/or form a closed curve. This is done by intro- 
ducing a parametric variable representation of the cubic spline function. Only cubic 
spline functions are dealt with in this report since they are naturally suitable to many 
engineering applications and offer a reasonable compromise between a simple and com- 
plex mathematical model. 


SYMBOLS 


Aj(x),Bj(x), 

Cj(x),Dj(x) 


coefficient functions for cubic polynomials 


c^,C 2 constants of integration 

E set of functions of coefficients 

F matrix of coefficients 

f(x) polynomial function in x 

G matrix of functions of coefficients 

g(t) polynomial function in t 

g(x) function describing equality of first derivative of adjoining cubics at jth joint 

h(t) polynomial function in t 

i index on elements in data sets 

j index of joints connecting cubic polynomials 

m number of joints 
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n number of input data points 

P degree of polynomial spline function 

R vector of residuals 

r residual 

S ordered data set in two dimensions 

S^(x) polynomial spline function of x over mesh A 

S A v (t) polynomial spline function of t over mesh A x 

Sa (t) polynomial spline function of t over mesh A v 

t independent variable for parametric technique 

t^ parametric variable function of x^ and y^ 

tj abscissa of jth element in set of joints for parametric independent variable 

tj abscissa of jth element in set of joints for parametric dependent variable 

W diagonal matrix of weights 

w weight 

X vector of ordinates of junction points and their second derivatives 

x independent variable for explicit technique or dependent variable for 

parametric technique 

Xi ith abscissa of data set S 

Xj abscissa of jth element in set of joints 

Y vector of ordinates of data points 
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y 


dependent variable for explicit technique 


ith ordinate in data set S 
y- ordinate of jth element in set of joints 

A set of interior mesh abscissas 

Ax, Ay sets of mesh abscissas for function of t 

A vector of Lagrange multipliers 

Xj Lagrange multiplier associated with jth spline 

cr standard deviation for explicit variable solution 

cr x ,Oy standard deviations for implicit variable solution 

<p function to be minimized 

7p constrained function to be minimized 

A prime indicates first derivative with respect to the independent variable; a double 
prime indicates second derivative with respect to the independent variable. 

DERIVATION 


The theory of least squares is applied to cubic polynomial spline functions S^(x) 
to derive the least -squares cubic polynomial fit for data smoothing. The least -squares 
development in a vector notation is reviewed in appendix A and the development of poly- 
nomial spline functions is discussed in detail in references 4 and 5. This report presents 
a brief development of the cubic polynomial spline function and the application of least 
squares to cubic splines. The explicit cubic polynomial splines are derived first, followed 
by a description of parametric splines and their relationship to explicit variables. 

Cubic Polynomial Splines 

A polynomial spline function is defined over an interval Xj ix = x n which is sub- 
divided into a mesh 
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A : Xj = xj < X2 < • • 


x j * * • <x m = x n 


Associated with this mesh is a set . . A polynomial spline function S^(x) is a set 

of m-1 polynomials of degree P connecting the points {(Xj,yj)} m ^ (called joints) 

such that the adjoining polynomials and the first P - 1 derivatives are continuous at the 
joints. 

A cubic spline function S^(x) is a set of cubic polynomials of degree three con- 
necting the mesh points {(x-j,y-j)} such that at the mesh points adjoining cubics and their 
first and second derivatives are continuous. This implies that S^(x), S^(x), and S^(x) 

are continuous in the interval xj £ x I x m . The derivation of a cubic spline fit to the 
set {(xj,y j )}! n i can be performed by first examining a single cubic in the interval 

xj ^ x = Xj + i and enforcing the continuity conditions with the adjoining cubic in the interval 
xj_l 2 x = xj. The following is a detailed explanation of the procedure. 

Let (xj,yj^ and (Xj + j,yj + i) be two joints (see fig. 1(a)) between which it is desirable 
to construct a cubic polynomial y(x). Since y(x) is a cubic polynomial, y"(x) (second 
derivative of y(x)) is a linear polynomial. The point slope representation (fig. 1(b)) for 
the second derivative is: 


y"(x) = 


y"( x j+i) 


y”(xj) 


c j+l 


Xi 


(x ■ x j) + y"( x j’ 


(i) 


Letting hj = Xj +1 - Xj, yj' = y"(Xj), and yj' +1 = y"(Xj + i), equation (1) can be rewritten 


y W - y j+l h-i 


+ y 3 


h i 


( 2 ) 


Integrating equation (2) twice, the following equations are obtained: 


„ x ( x “ x j) 2 ( x j+l " x ) 2 

y (x) = y j-n -lh f-- y i " - 2 ~ hj +C 1 


(3) 


( X ~ X 3 ) ■ ( X j+1 - x ) 

'3 6h{ 


yw - yj’+i - m t- * n ' +c i* +c i 


(4) 


3 
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(b) Second derivative. 

Figure 1. - Adjoining cubic polynomials 





Evaluating equation (4) at the jth and (j + l)th joints yields equations for cj and C 2 
as follows: 


yj+i 

c ! = — 







c 2-yj 


yj'hj- (yj'+i - yj>j*i (f j+i - fjft 
s 6 h i 


Collecting terms as coefficients of yj, yj', y.j + j, and yj' + j results in: 
y'(x) = Aj(x)yj + Bj(x)yj’ + Cj(x)yj +1 + D](x)yj’ +1 


(5) 


y(x) = Aj(x)fj + Bj(x)y-’ + Cj(x)yj +1 + Dj(x)yj’ +1 


(6) 


where 
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For each cubic there is one equation of the form of equation (6) and for m mesh 
points there will be m - 1 cubics. Having written the eubics in the form of equation (6), 
adjoining cubics at the interior joints %, . . ., x m _i} are equal and the second 

derivatives at the interior joints are equal. That is, 


y(Xj -) = y(xj +) 

y"(*j -) ■ y”(*j +) 


(j = 2, . . ., m - 1) 
(j = 2, , . ., m - 1) 


Enforcing the condition at the joints 

y’(*j -) - + ) 


yields the continuity of the spline function S^(x). Using equation (5) with j - 1 and j 
as indices establishes the above condition in the form 



Aj-l(xj)yj_i + [Cj-i(xj) - Aj.(xj|yj +Bj_ 1 (x. 


+ 


= 0 






(7) 


Since there are m - 2 interior points, there are m - 2 condition equations of the form 
of equation (7) that must be satisfied in order to have a cubic spline polynomial fit S^(x). 


Least -Squares Spline Solution 

The technique of least-squares polynomial curve fitting for smoothing noisy or 
irregular data represented by the data set S = ^x^'y^Vj 1 is standard. It consists of 
choosing a polynomial function and finding the polynomial coefficients which minimize the 
sum of the squares of the differences between the function values at the abscissas ^Xj}._j 
of the data set and the ordinates {yA n _ 1 of the data set. In the above, a set of cubic 
polynomials of the form of equation (6) with coefficients {(yj,yj')/j_j is defined on the 
interval jx^,x^| where the interval is divided into the mesh 

A : xj = xj < X2 < . , . xj . . . < x m - x n 
1 

In addition, at each ” the condition described by equation (7) must be satisfied. 

The resulting set of polynomials is a cubic spline function. 
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The following is a vector development of a least-squares solution procedure for 
determining the set {yj,y]'}j^l- Given a data set {( x i>yi)}£_i in the interval 
xi = x^xg x m = x n , the least-squares solution using equation (6) determines the 
set which minimizes 

f = Z r ? - f, [ y i - (8) 

i=l i=l 

Using the method of Lagrange multipliers and adjoining the conditions of equation (7), the 
constrained least-squares solution minimizes 


(p 


- J [n 

1=1 


n2 


m-1 


y ( x i) 




5=2 


( 9 ) 


The derivation of the condition that minimizes 7p is found in appendix A. This 
condition is expressed in matrix notation where the sets {(y^f an d 

((Yi - y(x i ))}"_ 1 are ordered into the column vectors: 


V 


’*1 

^2 


y i 

• 


y2 

* 

X = 

*% 

yi 


• 

• I 


^m 

yn 


fm_ 


yi - y( x i) 
y 2 - y( x 2 ) 


y n - y( x n) 


n ^ ^ n 

The sets {y j and {g( x j)}- =2 found by evaluating equations (6) and (7) at { x i} i=1 

and { x j}™2 are ex P ressed as FX and GX, respectively, where 
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Fj is a rectangular matrix with one row for each point in the interval xj lx 1 Xj + j., 
Sj is the number of points in the interval, and F is the composite of all the Fj's. 
There is coupling between Fj_ j, Fj, and Fj + j as indicated by the overlapping rec- 
tangles in F. Also, 



e 2,3 

e 2,4 

e 2,5 

e 2,6 



e 3, 1 

e 3,2 

e 3,3 

e 3,4 

e 3,5 

e 3,6 


Em-1,.1 Em-1,2 E m-1,3 E m _ 1>4 E m _ 1>5 
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and 
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Often, when working with a set of data, it is desirable to place greater emphasis on cer- 
tain points over others. This is achieved by associating a set of weights with 

the set such that each data point (x-i,y{) has a weight w^ and minimizing 


n m-1 

? - I ( r i w 0 2 + I x J g (*i) 

1=1 j=2 



where 



The primary mathematical tool in the form of equation (11) is now available. In the 
next section a step-by-step procedure to apply equation (11) to obtain satisfactory smooth 
fits to the set {( x i>yi)}^_i is presented. 

A Computational Algorithm for Explicit Variables y = f(x) 

Given the sets {(x^yi)}^, xj = xj and x m = x n , and {wj}^: 

(1) Choose 

(2) Determine the matrices F and G by computing the coefficients Ajj’Xj), 
Bj(Xi), Cj(xj), Dj(xi), and Ej ? k for k = 1, . . .,6. 
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(3) Form the matrices 






(6) Plot and the spline function S^(x). 

(7) If the resulting fit of S^(x) to the data is not satisfactory (i.e., 

<7 is not sufficiently small or S^(x) is not satisfactory from an engineering point of 
view), choose a new set and restart at (2). The number of joints m is not 

necessarily the same as in the previous attempt; however, the endpoints x^ and x m 
must still be equal, respectively, to and x n . Repeat until one of the following con- 
clusions is reached: 

(a) The fit is satisfactory and S^(x) can be used to represent 

(b) The fit is unsatisfactory and the technique is either abandoned or the 
parametric variable technique as described in the following section is applied. 

This approach is highly feasible in an interactive computer graphic environment. 
Rapid successive solutions can be attained and plotted for use in making a decision on the 
acceptability of a solution. Appendix B describes a computer program (D3670) written 
for the Control Data series 6000 computer systems using the CDC 250 series CRT display 
system and the LRC 6000 series graphics system. The technique can be applied in a non- 
interactive mode using off-line plotting equipment. 


Parametric Variables 

Using explicit variables (xj,y^) where y^ = f(x^) and finding the spline fit y = S^(x) 
as described in the previous sections required that on the mesh A, xj + i > xj. Con- 
sequently the described computer algorithm will not smooth arbitrary sets of data 
{( x i>yi)}^_j such as that shown in figure 2. If, however, a monotonic parametric 
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y 



Figure 2 . - Sample of data suitable for parametric fit. 

variable t is introduced to the data set such that the set becomes {(t^ xj), (t^yi)} ^ 
where = g(t^) and yj = h(ti), and choosing two meshes 


A x : ti = ti < t£ < . . • < t mi * t n 


A y : t l - ti < t 2 < . . . < tjjjg - t n 


the smoothing technique can be applied for arbitrary sets of data. A parametric variable 
which will satisfy the condition tj + i > tj while xj + j = xj is 

r o of/2 


r 2 2i 1/z 

*1+1 = |( x i+l - x i) + (n+1 - 7 i) J + *i 


ti = 0 


The accumulative chord length t will be monotonically increasing with respect to x 
and y provided (24+1,71+1) 4 - (x^yi). With the t's determined, there are two sets 
{( t b x i)}^ = l and {(*i>yi)}i_i upon which the algorithm described above can be applied. 
Applying this computational algorithm to the two sets yields the cubic spline functions 
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S Ax (t) and S Ay (t). For the purpose of interpolation, y can be found as a function of x 
by finding the inverse relation t = S Ax (x) and computing y = S Ay |$ Ax (x)j . The practical 

approach is to find t such that x - S Ax (t) = 0 in the interval jtj,tj + ij and substituting 
this t into y = S Ay (t). Since the cubic splines are being fitted with respect to a para- 
metric variable, the first and second derivatives of y as a function of x are: 

dy _ dy/dt _ S Ay^ 
dx dx/dt S Ax (t) 



For cases with data forming a closed curve, it is also desirable to have the spline 
curves continuous and equal in value at the endpoints. This is done by adding two addi- 
tional constraints in equation (9), thus introducing two additional Lagrange multipliers in 
the least-squares solution. Enforcing the condition at the endpoints 

y’(xi +) = y’(x m -) 

yields the continuity of the first derivative. Using equation (5) with 1 and m as indi- 
ces establishes the above condition in the form 

A i(*i)yi + B i(%)yi + c i(*i)y 2 +D i(*i)y2 

" V-l(%)yra-l * B m-l(%)ym-l “ C m _i(x m )y m - Dm-l(xm)ym = 0 

Enforcing equality at the endpoints requires the additional constraint 

fi - y m = 0 


A Computational Algorithm for Arbitrary Sets of Data 
Given the sets {[ x i,yi)}?_j and ^Wj)^: 

(1) Compute the set {ti}^_^ where tj = 0 and 

r o !/2 


r o 

ti+1 = [_(x i+ i - Xif + (y i+ i - yi ) j + ti 
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(2) Form the sets {(t^Xi)} ” =1 and {(t^y^^j. 


m - 1 

(3) Choose a set A x = {tj}._ ^ where tj = ti = 0 and t m ^ = t n . 


(4) Form the matrices 


(5) Find the inverse of 


'3 

F T WF g t_ 

. G 0 


and 


F t WF 


L G 


f t wf 


G J 


-lr 






F T WY| 


and solve 




f t wf 


(6) Compute standard derivation in the (t,x) space 

-n 1/2 


[ x i - 


n - 2m j 

V. 

J 


(7) Plot {(t i ,x i )}J 1 =1 and S Ax (t). 

(8) If the resulting fit of S Ax (t) to the data {(ti,x^ n _^ is not satisfactory, choose 

„ mi-1 

a new set (t-j). and restart at (4). Otherwise, continue. 

J 3=2 


(9) Choose a set Ay = {tj} ^ where tj = t^ and t m ^ = t n . 

a f| 


-1 


(10) Form the matrices 


(11) Find the inverse of 


F T WF G 

G 0 

F T WF 


and 


J t 


f t wy 


and solve 


A 

■- J y 


f t wf g t 

G 0 


J t 


f t wf 
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(12) Compute standard derivation in the (t,y) space 

1/2 



n - S Ay(*i)] 

2 w?l 

n - 2 ni 2 



(13) Plot {(tpYi)}^ and S Ay (t). 

(14) If the resulting fit of S A (t) to the data {(tpYi)} 11 , is not satisfactory, choose 

m 1 ^ • 

a new set {ti} m ^ and restart at (10). Otherwise, continue. 
iJ ]=2 

(15) Plot {(Xi, yi )}“ =1 and S Ay (S Ax (t)). 

(16) If the resulting fit of S Ay ^S Ax (t)j to is not satisfactory, restart 

at (3). Repeat until one of the following conclusions is reached: 

(a) The fit is satisfactory so that S Ax (t) and S Ay (t) can be used to repre- 
sent {(x i ,y i )}^ r 

(b) The fit is unsatisfactory and the technique is abandoned. 

Here again the technique is highly feasible in an interactive computer graphics envi- 
ronment. In appendix B, the described computer program has the option of using para- 
metric variables in the interactive and off-line modes. In appendix C there is a descrip- 
tion of example cases using both explicit and implicit variables. 


CONCLUDING REMARKS 


Two algorithms have been developed for smoothing sets of data in two-dimensional 
Cartesian coordinates. Both algorithms exploit the cubic spline function which is contin- 
uous and has continuous first and second derivatives. A least-squares technique is applied 
to the cubic spline function to obtain smooth representations of the data sets. The two 
algorithms, the explicit variable algorithm and the parametric variable algorithm, are 
distinguished according to how the independent and dependent variables are defined rela- 
tive to each other. The explicit variable algorithm requires less computation and is, 
therefore, recommended whenever high gradients are not encountered. The parametric 
variable algorithm requires more computation and can be applied to arbitrary sets of data. 
Both algorithms require that the user specify a set of mesh points in the independent 
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variable direction. Since this mesh usually must be modified in order to achieve satis- 
factory solutions, the algorithms are best applied in an interactive computer graphics 
environment. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., November 1, 1973. 
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APPENDIX A 


CONSTRAINED WEIGHTED LINEAR LEAST -SQUARES ESTIMATION 
Given the data set an( j the function 

y = y(x, y 1} y{, , . y m , y^) ( 2m < n ) 

a weighted least -squares solution is the set of estimates 

(«2k-i = yk> a 2k = yk}“ =1 

which minimizes the sum of the squares of the weighted differences between the corre- 
sponding elements of the sets ^y^jV^ and {y(xt, ap a 2 , . * ( see re E ®).» 

that is, 



If there exists a set of functions 

a 2> • • •> a i)}p-i 


(l + q < n) 


such that 

*“<%>(«!» a 2> • • •> a l ) * °}p = i 

then if/ defines a set of constraints which the function y(x, ap a 2, . ■ a^) must satisfy. 

The constrained least-squares problem consists of determining the set {®k}k = i which 
minimizes equation (Al) subject to the conditions described by if/. This is accomplished 
by associating a set of Lagrange multipliers /\nV* with the set \f/ and determining 
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APPENDIX A - Continued 



An additional mathematical assumption is that 
y(x, a h a 2 , . . a t ) 

and 

{gp(«P a 2> • • •> a i)}p =1 

are linear functions with respect to the elements of the set an d, therefore, can 

be written in the form 

l 

y(x, a h a 2 , , . <xj) = J a k (x)o! k 

k=l 

l 

gp( a l> a 2> • • •* a l) = X b P kQ!k 

k=l 


With these definitions, the minimization problem for constrained weighted linear least - 
squares estimation can be rewritten in matrix notation as 

min $ = min|R T WR + A T BX 
X, A X, A 

= min (y T - X T A T )w(Y - AX) + 2A T BX (A3) 

X,A L j 

where 


n 


r l 


a l( x l) 

a 2( x l) 

• • • a z( x i> 

y2 

R = 

r 2 

A = 

a l( x 2) 

a 2( x 2) 

• • • a z( x 2) 

y n 


r n_ 


a l( x n) 

a 2( x n) 

• • • a i( x n) 
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APPENDIX A - Continued 


«1 b ll b 12 • * • b U H 

a 2 b 2 l b 22 • • • b 2 1 x 2 

X= • B = ’ ' A = ’ 

• • • • • 
oil bqi bq2 . • . bq^ Xq 



To determine min the first variation 
X, A 

that is, 

5 $ = ~2(y T WA - X T A T WA) 
Since 6X * 0 and 6 A * 0, 

-A T WY + A T WAX + B T A = 0 
BX = 0 



6$ with respect to X and A must vanish; 
A T B5X + 26A T BX =0 

(A4a) 

(A4b) 


This can be written in a single matrix equation as 



(A5) 
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or 

A = b[a T Wa] 

T* 

The vector A has a unique solution if A WA is invertible and if B 
invertible. If these conditions are satisfied for A, then X has a unique solution by sub- 
stituting A into equation (A7). 

In addition, a sufficient condition that solution of equations (A4) yields a unique global 
minimum of equation (A2) subject to the constraints BX = 0 is that SX^A^WASX > 0 for 
all 5X > 0 with B6X = 0 (ref. 7). 
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APPENDIX B 


PROGRAM DESCRIPTION 
General Discussion 

The program, written in FORTRAN IV for the CDC 6000 series computers, consist^ 
of a main program and five subroutines. CalComp plotting routines and CDC 250 routines 
are used for CRT display. The program optionally can run in either the batch mode or on 
line with control from a CRT console. 

The main program directs the solution procedure, either from preset instructions in 
the batch mode or according to commands from the CRT console in the on-line mode. In 
the on-line mode, instructions to the user are displayed to indicate what options are avail- 
able at the various points in the program. 

The following subroutines are called from the main program or its subprograms. 

Subroutine Function 

PLPT scales and plots computed values with the CalComp POINT routine 

CUSPFIT determines a cubic equation approximating the input data and computes the 
first and second derivative 

SUP computes the parametric variable t and sets up the arrays for subroutine 

CUSPFIT in cases using the parametric option of the program 

MINMAX finds the minimum and maximum values of the data to be plotted 

SCALEBW scales the values for plotting in the parametric cases 

SIMEQ solves the matrix equation AX = B where A is a square coefficient 

matrix and B is a matrix of constant vectors 

FTLUP performs a second-order interpolation to find intermediate values from a 

tabular array 

Interactive Graphics and Plotting Routines 

The CRT routines CDC 250, SCREEN, PARAMS, MESSAGE, and NEXT are from the 
LRC interactive graphics software package and the plotting routines CALCOMP, LEROY, 
ASCALE, AXES, LINPLT, LINE, NOTATE, NUMBER, CALPLT, and INFOPLT are from 
the CalComp software package. Plotter output is routed to a tape during job execution and 
after job completion is plotted on a CalComp digital incremental plotter. The arrays X, 
XC, Y, YY, YJOIN, COMPY, and R are used for plotting. 
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APPENDIX B - Continued 


Descriptions, Flow Charts, and Listings of the Main 
Program and Subprograms 

Main program . - The main program directs the solution procedure from preset 
instructions in the batch mode, or from CRT console commands in the on-line mode. 
The flow diagram of the main program is as follows: 
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APPENDIX B - Continued 


The program listing for the main program is as follows: 


* 

* 

* 

* 

* 

* 

* 

* 


* 

* 

♦ 


4 


1 


2 


3 

4 


5 


6 


* 


* 


PROGRAM MA INC INPUT =201 * OUT PUT = 201 ,T APE 5= INPUT* T AP£6=0DTPUT) 

A 

I 

100000 

200000 

X AND YY ARE COORDINATES OF THE INPUT DATA POINTS 

A 

2 

300000 

R-LEFT ENDPOINT, JUNCTION PO INTS* RIGHT ENDPOINT 

A 

3 

400000 

NK RENUMBER OF CURVES 

A 

4 

500000 

NMAT- NUMBER OF POINTS 

A 

5 

600000 

IFLAG =0 DATA POINT COORDINATES ARE READ WITH FORMAT 

A 

6 

700000 

IFLAG =1 DATA POINT COORDINATES ARE READ WITH NAMELIST 

A 

7 

800000 

w - ARRAY of weights of data points 

A 

8 

900000 

KO = 1 STANDARD CRT 

A 

9 

1000000 

KO = 2 STANDARD BATCH 

A 

10 

1100000 

KO = 3 PARAMETRIC CRT 

A 

11 

1200000 

KO = 4 PARAMETRIC BATCH 

A 

12 

1300000 

KLOSF * o ARBITRARY CURVE 

A 

13 

1400000 

KLOSE = I CLOSED CURVE (MUST BE PARAMETRIC) 

A 

14 

1500000 

1600000 

REAL MAXRES 

A 

15 

1700000 

DIMENSION TSAVEC3I0) 

A 

16 

1800000 

DIMENSION CS { 54*4) * Z(310), 7(31(1), RS<30)* 0130) 

A 

17 

1900000 

DIMENSION XS AV El 310 i • YSAVE(3iO> 

A 

18 

2000000 

COMMON XRAY< 3) ,YRAY<3> 

A 

19 

2100000 

COMMON X1310),Y(310), A( 304,4) ,R( 30) , C ( 84 . 85 1 ,L< 841 »N( 28) , YYI 3 1.01 

A 

20 

2200000 

COMMON COMPY ( 310 ) » W( 3 10 ) , Y JOIN! 30 ) , XC ( 310 ) 

A 

21 

2300000 

NAMELIST /NAMI/ X, YY ,NKR , R ,NM AT , I FL AG . W, KO , KLOSE 

A 

22 

2400000 

2500000 

INITIALIZE PLOT ROUTINE AND SET ORIGIN 

A 

23 

2600000 

2700000 

DATA NKR/i/, IFLAG/ 1/ , KO/ 1/ • NMAX /304/, KLOSE/O/ 

A 

24 

2800000 

IBEN=0 

A 

25 

2900000 

CONTINUE 

A 

26 

3000000 

R( U =0. 

A 

27 

3100000 

DO 2 1=1,310 

A 

28 

3200000 

w(ii=i. 

A 

29 

3300000 

IL = l 

A 

30 

3400000 

wr=i*o 

A 

31 

3500000 

XPG=10. 

A 

32 

3600000 

YPG= l 0 • 

A 

33 

3700000 

X0V=10. 

A 

34 

3800000 

NC=84 

A 

35 

3900000 

XT IC=1 • 

A 

36 

4000000 

Y Q V= 1 0 * 

A 

37 

4100000 

YT I C= l « 

A 

38 

4200000 

IF CIBEN.EQ.O) GO TO 4 

A 

39 

4300000 

DO 3 1=1, N MAT 

A 

40 

4400000 

XI I ) =X SAVE < I ) 

A 

41 

4500000 

CO NT INUE 

A 

42 

4600000 

READ ( 5 , NA Ml ) 

A 

43 

4700000 

IF I IBEN*EO* 1) GO TO 7 

A 

44 

4800000 

IBEN= I 

A 

45 

4900000 

GO TO (6*5, 6,5) , KO 

A 

46 

5000000 

CALL CALCOMP 

A 

47 

5100000 

CALL LEROY 

A 

48 

5200000 

GO TO 7 

A 

49 

5300000 

CALL CDC250 

A 

50 

5400000 

5500000 

SETS PARAMETERS IN PLOT ROUTINE TO OUTPUT OR CHANGE AT CRT 

A 

51 

5600000 

5700000 

CALL SCREEN II., 1., .9) 

A 

52 

5800000 

CALL PARAMS 

A 

53 

5900000 

6000000 

THIS CLEARS THE PARAMETER TABLE 

A 

54 

6100000 

6200000 

CALL PARAMS I2LAN, AN, 3LSTD,ST0) 

A 

55 

6300000 

CALL PARAMS (2LRI, R< 1) * 2LR2,R ( 2 ) , 2LR3,R( 3) ) 

A 

56 

6400000 

CALL PARAMS (2LR4, R(4) ,2LR5.R ( 5 ) ♦ 2LR6, R ( 6) ) 

A 

57 

6500000 

CALL PARAMS ( 2L R7, R < 7> , 2LR 8, R( 8 ) ,2LR9 ,Rl 9) ) 

A 

58 

6600000 

CALL PARAMS ( 3LRiO,RMO) ,3LR11 , Rf 11 ) > 

A 

59 

6700000 
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APPENDIX B - Continued 



CALL P ARAMS 1 3LR12,R{ 12) t3LR13, R( 13) » 3LR14*RL14) ) 

A 

60 

6800000 


CALL P ARAMS i 3LR15 ♦ R ( 15) ,3LR16, R( 16 ) , 3LR17.Rf 17 )) 

A 

61 

6900000 


CALL PARAMS { 3LR18,RC 18) , 3LR19,R( 191 ,3LR2Q ,R(20)) 

A 

62 

7000000 


CALL PARAMS i 3LR21 , R{ 21 ) ) 



A 

63 

7100000 


CALL PARAMS < 2LA l, Al, 2LWT , WT , 1L l , U 



A 

64 

7200000 


CALL MESAGE U,39HPR0GRAM U329Q —FINDS 

BEST 

FIT FOR DATA, 391 

A 

65 

7300000 


CALL MESAGE (1,32HUSE BEST FIT TO SOLVE 

PARAMETERS, 321 

A 

66 

7400000 


CALL MESAGE (1,40HHIT KEY 46 TO READ DATA TO 

START PROGRAM, 401 

A 

67 

7500000 


CALL MESAGE <1,22HHLT KEY 45 TO STOP JOB 

,22) 


A 

68 

7600000 


CALL NEXT iNKEYJ 



A 

69 

7700000 


IF INKEY.NE.46I GO TO 48 



A 

70 

7800000 

7 

IF (KQ.LT.3) K0DE=4 



A 

71 

7900000 


IF (K0.E0.2I NKEY=49 



A 

72 

8000000 


IF LK0.EQ.4) NKEY=51 



A 

73 

8100000 


IF UFLAG.EQ.O) GO TO 9 



A 

74 

8200000 


00 8 1=1* NMAT 



A 

75 

8300000 


YS AVE I n = YY( I) 



A 

76 

8400000 

8 

xsavei i)=xm 



A 

77 

8500000 

9 

IF ( KO.GT .21 KODE= 0 



A 

78 

8600000 


AN=NKR 



A 

79 

8700000 


NKRi=NKR+i 



A 

80 

8800000 


IF (ENDFILE 5) 47,10 



A 

81 

8900000 

13 

IF (IFLAG.NE.O) GO TO 16 



A 

82 

9000000 


NMAT=G 



A 

83 

9100000 

11 

NMAT =NMAT+ 1 



A 

84 

9200000 


READ (5,75) X( NMAT ), Y( NMAT ) 



A 

85 

9300000 


IF ( X ( NMAT >.NE. lliill.) GO TO 11 



A 

86 

9400000 


GO TO 113, 13,12,12), KO 



A 

87 

9500000 

12 

CONTINUE 



A 

88 

9600000 


CALL SUP < X.Y,T,Z,XMAXT,NMATI 



A 

89 

9700000 


XM=9H T VALUES 



A 

90 

9800000 


YM = 9H X VALUES 



A 

91 

990.0000 


GO TO 14 



A 

92 

10000000 

13 

XM=9H X VALUES 



A 

93 

to 100000 


YM=9H Y VALUES 



A 

94 

10200000 

14 

CONTINUE 



A 

95 

10300000 


NMAT=NMAT-l 



A 

96 

10400000 


DO 15 1=1, NMAT 



A 

97 

10500000 

15 

YY1 1 )=Y Cl) 



A 

98 

10600000 


GO TO 21 



A 

99 

10700000 

16 

GO TO (19,19,17,171, KO 



A 

100 

10800000 

17 

CALL SUP ( X, YY »T , Z , XMAXT , NMAT i 



A 

101 

10900000 


XM-9H T VALUES 



A 

102 

11000000 


YM=9H X VALUES 



A 

103 

11100000 


DO 18 1=1, NMAT 



A 

104 

11200000 

18 

TSAVEi I) = Xt I) 



A 

105 

11300000 

19 

DO 20 1=1, NMAT 



A 

106 

11400000 

20 

Y( 1 1 = Y Y( I'l 



A 

107 

11500000 


IF (K0.GE.3) GO TO 21 



A 

108 

11600000 


XM=9H X VALUES 



A 

109 

11700000 


YM=9H Y VALUES 



A 

110 

11800000 

21 

NKR1=N KR+1 



A 

111 

11900000 


IF lKQ.jBE.3l GO TO 23 



A 

112 

1200*0000 


XO=X( 1 ) 



A 

113 

12100000 


DO 22 1=1, NMAT 



A 

114 

12200000 

22 

xm=xm-xo 



A 

115 

12300000 


Rm=xcn 



A 

116 

1*2400000 

23 

CONTINUE 



A 

117 

12500000 


IF IKODE.EO.OI WRITE 16,741 



A 

118 

12600000 


IF (KOOE.EQ.i) WRITE 16,81) 



A 

119 

12700000 


IF (K00E.EQ.4) WRITE (6,88) 



A 

120 

12800000 


RCNKR1 )=X( NMAT) 



A 

121 

12900000 


WRITE (6,76) NKR, NMAT , (R ( J ) , 1= 1 , NKR l ) 



A 

122 

13000000 


LL=NKR*2+2 



A 

123 

13100000 


LR=LL+NKR 



A 

124 

13200000 


LA=LR- 1 



A 

125 

13300000 


LO=LA 



A 

126 

13400000 


CALL CUSPFIT U , C, NKR , L, X, Y, R , NMAX ,NC ,N, 

> A, W, 

KODE.CS.KLOSE) 

A 

127 

13500000 


I S TUNT=N( NKR + 1 ) 



A 

128 

13600000 
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CALL CUSPFIT <2,C,NKR,L,X,Y,R, NMAX,NC*N* A, W,KQQE,CS, KLQSEI A 129 

RES=0 A 130 

IF { KODE* E 0* 0) WRITE (6,78) A 131 

IF { KQDE.EQ. I ) WRITE (6,77) A 132 

IF (KQDE.EQ. 4) WRITE (6,87) A 133 

DO 24 1=1, ISTUNT A 134 

CGMPYl D=Y( I) A 135 

RES I D= YY 1 1 )- Y( 1 ) A 136 

IF (W(l).EO.O.) RE S I D= 0 • A 137 

WRITE (6,791 t, XI i )# YY( I ) ,COHP Y( l )» RE SID A 138 

RE5=RESID**2*W< II A 139 

24 CONTINUE A 140 

RES=RES/( I STUNT-2*NKRI A 141 

ST D= SORT (RES) A 142 

PRINT 80, STD A 143 

N0= I STUNT A 144 

[F (K0.GT.2) GO TO 28 A 145 

* PUT IN MORE POINTS IF ORIGINALLY LESS THAN 50 FOR STANDARD VERSION A 146 

IF (N0.GT.5Q) GO TO 27 A 147 

NM 1=NQ— 1 A 148 

NP=120/NM1 A 149 

NPM1=NP-1 A 150 

I C=1 A 151 

XC(l)=X(l) A 152 

DO 26 l 1=1 ,NM1 A 153 

FACs(X(I.lFU-X( im/NPMl A 154 

DO 25 [=1 , NPM1 A 155 

IC=IC*1 A 156 

-XCUC)=xC.! IC-ll+FAC A 157 

25 CONTINUE A 158 

26 CONTINUE A 159 

XCUC) =R(NKRU A 160 

CALL CUSPFIT l 2 , C. NKR, L , XC , COMPY, R, NMAX, NC , N, A , W ,KODE , CS , KLO SE I A 161 

27 CONTINUE A 162 

28 CALL CUSPF IT I 5 , C * NKR ,L , X , Y,R , NHAX » NO * N* A, W, KODE, CS • KLOS E ) A 163 

K= 1 A 164 

FNKR=NKR A 165 

* COMPUTE MINI MUMS AND MAXIMUMS A 166 

CALL ASCALE ( X , XPG . NO, K , 10 - ) A 167 

CALL ASCALE ( YY , YPG, NO , K , 1 0. ) A 168 

* .DRAM X AXIS A 169 

NPi=NO+l A 170 

NP 2= NO + 2 A 171 

CALL AXES 10. , 0 . , 0 . . XPG, X I NP 1 ) • X( NP 2) , XT 1C , XDV , XM, . 15,-9) A 172 

* DRAW Y AXIS A 173 

CALL AXES (0 •* 0. ,90. » YPG, YY( NJP IT * YY(NP2) , YTIO, YDV, YM» • 15,91 A 174 

* PLOT CURVE A 175 

CALL PLPT (X, YY, NO) A 176 

NP=I C A 177 

IF (K0*LT,3.AND.N0.LT.50) GO TO 30 A 178 

NP=NQ * A 179 

DO 29 1=1, NO A 180 

29 XC( I ) = XU) A 181 

30 NN 1=NP 4-1 A 182 

NN2=NP+2 A 183 

XCI NN1 )=X ( NP 1) A 184 

XC ( NN2 )=X( NP2) A 185 

COMPY ( NN1 )=Y Y( NP 1 ) A 186 

COMP Y( NN2 ) =YY ( NP2) A 187 


13700000 
13800000 
13900000 
14000000 
14100000 
14200000 
14300000 
14400000 
14500000 
14600000 
14700000 
14800000 
14900000 
15000000 
15100000 
15200000 
15300000 
15400000 
15500000 
15600000 
15700000 
15800000 
15900000 
16000000 
16100000 
16200000 
16300000 
16400000 
16500000 
16600000 
16700000 
16800000 
16900000 
17000000 
17100000 
17200000 
17300000 
17400000 
17500000 
17600000 
1770 0000 
17800000 
17900000 
18000000 
18100000 
18200000 
18300000 
18400000 
18500000 
18600000 
18700000 
18800000 
18900000 
19000000 
19100000 
19200000 
19300000 
19400000 
19500000 
19600000 
19700000 
19800000 
19900000 
20000000 
20100000 
2020.0000 
20300000 
20400000 
20500000 


29 



APPENDIX B - Continued 


31 


32 


* 


33 


* 


34 


35 




36 


CALL LINPLT CXC, CO MPY*NP*K,0, 0*0,0) 

A 

188 

20600000 

NK l i=N KR 1+ 1 

A 

189 

20700000 

NK12=NKR 1 + 2 

A 

190 

20800000 

RI NK1 l ) — X< NP 1 i 

A 

191 

20900000 

R<NK12)=X<NP2) 

A 

192 

21000000 

YJOINC NK1 1 ) = YYC NPl 1 

A 

193 

21100000 

YJOI N{ NKL2 1 “YY INP2 1 

A 

194 

21200000 

YJQl N< l)=COMPY( i ) 

A 

195 

21300000 

DO 31 1=2 * NKR1 

A 

196 

21400000 

CALL FTLUP C R € 1 1 « YJOINC L), 1, NP* XC»CGHPY) 

A 

197 

21500000 

CALL L INE ( R , Y JO I N *NKR1 , Kt— 1 ,5**5) 

A 

198 

21600000 

CALL NOTATE U - , 9- 5, . 14* 19 HNJHBER OF CURVES » , 0*,19) 

A 

199 

21700000 

CALL NUMBER t 4* , 9. 5# . 14* FNKR, 0. ,-l) 

A 

200 

21800000 

CALL NOTATE <5 • • 9* 5* *14, 32HENQPGINT $ AND JUNCTION POINTS * ,0*,32) 

A 

201 

21900000 

XN = 0 • 

A 

202 

22000000 

YNN=9* 0 

A 

203 

22100000 

DO 32 l-l.NKRi 

A 

204 

22200000 

XN*XN+.7 

A 

205 

22300000 

IF 1 I.NE.14) GO TO 32 

A 

206 

22400000 

YNN=8.5 

A 

207 

22500000 

XN=*7 

A 

208 

22600000 

CALL NUMBER ( XN» YNN, • 1 1«R{ 11 • 0. *4) 

A 

209 

22700000 

22800000 

ESTABLISH A NEW REFERENCE POINT FOR THE NEXT GRAPH 

A 

210 

22900000 

23000000 

CALL CALPLT 10..0..-3I 

A 

211 

23100000 

GO TO 4 33,34,33, 34 ) , KD 

A 

212 

2 3200000 

CONTINUE 

A 

213 

23300000 

23400000 

THIS WILL STOP PROGRAM 

A 

214 

23500000 

23600000 

CONTINUE 

A 

215 

23700000 

CALL CALPLT U2.,0.,-3) 

A 

216 

23800000 

IF (K0.E0.2) GO TO 49 

A 

217 

23900000 

IF (IL.EQ.2) CALL BALLPT 

A 

218 

24000000 

IL = 1 

A 

219 

24100000 

IF <K0*EG*4) GO TO 41 

A 

220 

24200000 

CALL MESAGE ( 1 ,44HK4 5-CHANGE WEIGHTS K46-SET ALL WE IGHTS=1 *44) 

A 

221 

24300000 

CALL MESAGE C1.23HHIT ANY KEY TO CONTINUE, 23) 

A 

222 

24400000 

CALL NEXT (NKEY) 

A 

223 

24500000 

IF CNKEY.EQ.45) GO TO 36 

A 

224 

24600000 

IF (NKEY.NE.46) GO TO 39 

A 

225 

24700000 

DO 35 1=1, NO 

A 

226 

24800000 

WC IT =1.0 

A 

227 

24900000 

GO TO 39 

A 

228 

25000000 

25100000 

LOOK AT EACH POINT AND WEIGHT THEM IF YOU WANT TO 

A 

229 

25200000 

25300000 

CALL PLPT (X * YY, NO ) 

A 

230 

25400000 

XRA Y ( 2 ) = X I NPl) 

A 

231 

25500000 

XR AY ( 3 )=X ( NP2 ) 

A 

232 

25600000 

YRAY 4.2 i -YY lNPl J 

A 

233 

25700000 

YRAYI3MYY INP2) 

A 

234 

25800000 

All=0. 

A 

235 

25900000 

A1=0. 

A 

236 

26000000 

DO 37 11=1, NO 

A 

237 

26100000 

A11=A1 1+Al 

A 

238 

26200000 

1= II FA 11 

A 

239 

26300000 

A1 = 0* 

A 

240 

26400000 

XRAY ( 1 ) = XC I ) 

A 

241 

26500000 

YRAY ( 1 )=YY (I 1 

A 

242 

26600000 

CALL LINE i XRAY , YRAY , l ,K 1, 4, *09) 

A 

243 

26700000 

CALL CALPLT CO* , 0* ,-3J 

A 

244 

26800000 

CALL MESAGE (l,48HANY K CONTINUE K 4 5-CHANGE WT K46-ST0P SEARCH 

A 

245 

26900000 

1,48) 

A 

246 

27000000 

CALL MESAGE (l,35HWT = NEW WEIGHT A1 = SKIP POINTS, 35) 

A 

247 

27100000 

CALL NEXT INKEY) 

A 

248 

27200000 

IF C NKEY. E 0* 46 ) GO TO 38 

A 

249 

27300000 

IF (NKEY.NE.45) GO TO 37 

A 

250 

27400000 

27500000 
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♦ 


CHANGE WEIGHT OF POINT 

A 251 

27600000 





27700000 



W( I ) =WT 

A 252 

27800000 


37 

CONTINUE 

A 253 

27900000 


38 

CALL CALPLT U2. »0.*-3l 

A 254 

28000000 


39 

CONTINUE 

A 255 

28100000 



CALL MESAGE { i , 42HK45-T YPE NEW R K 46- NEW CASE K47-ST0P.42) 

A 256 

28200000 



CALL MESAGE (1.30HK-49 TO PLOT FIRST DER I VAT IVES .301 

A 257 

28300000 



IF IKO.NE.3I GO TO 40 

A 258 

28400000 



CALL MESAGE U# 19HK-44 TO PLOT Y VS ft 191 

A 259 

28500000 



CALL MESAGE U.19HK-50 TO PLOT Y VS X.19) 

A 260 

28600000 



CALL MESAGE U.34HK-51 TO READ NEW R AND PLOT Y VS T.34) 

A 261 

28700000 



CALL MESAGE (U25HK-52 TO RE-COMPUTE X VS T.25) 

A 262 

28800000 



CALL MESAGE U.25HK-53 TO RE-COMPUTE Y VS T.25) 

A 263 

28900000 





29000000 

* 


THIS STOPS PROGRAM 

A 264 

29100000 





29200000 


40 

CALL NEXT (NKEY) 

A 265 

29300000 


41 

NKR= AN 

A 266 

29400000 



IF CKODE.NE.il GO TO 43 

A 267 

29500000 



Ni=NKRl 

A 268 

29600000 



NK= NKR 

A 269 

29700000 



00 42 I- L t NKR1 

A 2 70 

2980.0000 


42 

RS Cl I — R 1 1 .) 

A 271 

29900000 


43 

CONTINUE 

A 272 

30000000 



XI I STUNT* 1 ) = R i NKR* l 1*1. 

A 273 

30100000 



IF (NKEY. EG. 45) GO TO 19 

A 274 

30200000 



IF C NKEY.E 0. 44 ) GO TO 50 

A 275 

30300000 



IF (NKEY.E0.50i GO TO 57 

A 276 

30400000 



IF (NKEY.E 0.511 GO TO 50 

A 277 

30500000 



IF (NKEY.E 0.461 GO TO 1 

A 2 78 

30600000 



IF INKEY.E0.52) GO TO 44 

A 2 79 

30700000 



IF (NKEY.E0.53i GO TO 50 

A 280 

30800000 



IF (NKEY.E 0.49) GO TO 49 

A 281 

30900000 



IF ( NKEY.E0.47 ) STOP 

A 282 

31000000 



GO TO 39 

A 283 

31100000 





31200000 

* 


RE-COMPUTE X VS. T 

A 284 

31300000 





31400000 


44 

DO 45 1 = 1 ♦ NMAT 

A 285 

31500000 



X( n^TSAVEl l i 

A 286 

31600000 



YY.C1 ) = XSAVEI l) 

A 287 

31700000 



T4 U =T SAVE ( I i 

A 288 

31800000 


45 

Z( I J =Y SAVE ( I 1 

A 289 

31900000 



XM=9H T VALUES 

A 290 

32000000 



YM=9H X VALUES 

A 291 

32100000 



IF (K0.GT.2) KODE=0 

A 292 

32200000 



DO 46 1=1* NKRS 

A 293 

32300000 


46 

R ( I ) =Q< l ) 

A 294 

32400000 



NK=NKR 

A 295 

32500000 



N1-NKR1 

A 296 

32600000 



NKR=NKS 

A 297 

32700000 



NKR1=NKRS 

A 298 

32800000 



GO TO 19 

A 299 

32900000 


47 

CALL CALPLT (0.0*999) 

A 300 

33000000 


48 

STOP I 

A 301 

33100000 





33200000 

* 


PLOT FIRST DERIVATIVES 

A 302 

33300000 





33400000 


49 

CALL CUSPF IT (3. C» NKR* L. XC . COMP Y* R* NMAX . NC *N*A. W.KODE, CS. KLOSEi 

A 303 

33500000 



CALL ASCALE JXC ,XPG, NP *K ,1 0. ) 

A 304 

33600000 



CALL ASCALE ( COM PY . YPG. NP *K . LO I 

A 305 

33700000 



CALL AXES (0..0. , 0. * XPG. XC i NN1 ) »XC( NN2 i * XT IC.XDV *XH* • 15* — 9i 

A 306 

3380.0000 



CALL AXES ( 0 • * 0. * 90. * YPG. COMPY ( NN 1) , COMP YC NN2) . YTIC * YD V, YM 15,9 • 

A 307 

33900000 



CALL LINPLT { XC . COMPY. NP ,K .0 ,0 , 0, 0) 

A 308 

34000000 



CALL CUSPF IT 1 3. C, NKR. L . X* COMPY, R. NMAX.NC, N* A.W , KODE ,C S.KLOSE ) 

A 309 

34100000 



WRITE (6.72) 

A 310 

34200000 



WRITE (6.731 ( X ( I ) . YY ( I ) .COMPY ( I ) » 1=1 . NO ) 

A 311 

34300000 
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CALL CALPLT U2**0.,-3i 

A 

312 

34400000 


IF (KQ.EQ.2) GO TO 1 

A 

313 

34500000 


GO TO 39 

A 

314 

34600000 

50 

DO 51 1-1 * NMAT 

A 

315 

34700000 

34800000 

* 

COMPUTE Y VS* T 

A 

316 

34900000 

35000000 

51 

YYl I ) = Z(M 

A 

317 

35100000 


YM=9H Y VALUES 

A 

318 

35200000 


XM=9H T VALUES 

A 

319 

35300000 


IF { KO..«-.GT *2) KQDE= 1 

A 

320 

35400000 


IF (NKEY.EQ.53) GO TO 54 

A 

321 

35500000 


DO 52 l~l,NKRl 

A 

322 

35600000 

52 

om=Rm 

A 

323 

35700000 


NKRS=NKR1 

A 

324 

35800000 


NKS=NKR 

A 

325 

35900000 


IF INKEY.EQ.5i) READ (5.NAM1) 

A 

326 

36000000 


IF (ENDFILE 5) 47,53 

A 

327 

36100000 

53 

AN=NKR 

A 

328 

36200000 


IF ( KO *FQ* 4) NKEY= 50 

A 

329 

36300000 


IF INKEY.NE.53) GO TO 19 

A 

330 

36400000 

54 

00 55 1=1, N1 

A 

331 

36500000 

55 

R( I )=RS( 11 

A 

332 

36600000 


00 56 1= 1* NMAT 

A 

333 

36700000 


XI n = TSAVE(I) 

A 

334 

36800000 

56 

YY ( I ) = YSAV E( I) 

A 

335 

36900000 


NKR=NK 

A 

336 

37000000 


NK R1=N 1 

A 

337 

37100000 


GO TO 19 

A 

338 

37200000 

37300000 

♦ 

COMPUTE X VS.Y 

A 

339 

37400000 

37500000 

57 

Tin=o. 

A 

340 

37600000 


XM=9H X VALUES 

A 

341 

37 7000 00 


DT=XMAXT / 300 • 

A 

342 

37800000 


DO 58 1=1 , Hi 

A 

343 

37900000 

58 

Rl 1 )=RSC I ) 

A 

344 

38000000 


NKR1=N1 

A 

345 

38100000 


NKR=NK 

A 

346 

38200000 


'DO 61 1=2*301 

A 

347 

38300000 


T ( I)=TI l-U+DT 

A 

348 

38400000 


00 59 J= 1 * NKR1 

A 

349 

38500000 


J J = NKR 1— J + 1 

A 

350 

38600000 


JI = JJ 

A 

351 

38700000 


IF (Tin.GE.R(JJ)l GO TO 60 

A 

352 

38800000 

59 

CONTINUE 

A 

353 

38900000 

60 

IF IJl.GT.NKR* Jl = NKR 

A 

354 

39000000 

61 

YC I) = ( 1 CSI J I + 27 » 1 ) *T ( I )*CS( JI + 2 7* 2) )*TII H-CSl JI #-27 »3) ) *T Cl )+CS( J1 + 

A 

355 

39100000 


127*4) 

A 

356 

39200000 


DO 64 1=2*301 

A 

357 

39300000 


DO 62 J=1 * NKRS 

A 

358 

39400000 


JJ=NKRS-J+1 

A 

359 

39500000 


JI = JJ 

A 

360 

39600000 


IF ITMI.GE.Ql JJ >) GO TO 63 

A 

361 

39700000 

62 

CONTINUE 

A 

362 

39800000 

63 

IF (JI.GT.NKS) J I=NKS 

A 

363 

39900000 


xc n=i < esc ji *d*t( n*csi ji ,211 *rm+csu 1* 31 t*tc n +csc ji*4) 

A 

364 

40000000 

64 

CONTINUE 

A 

365 

40100000 


xi n=csu*4) 

A 

366 

40200000 


Y( 1 )=CS (28*4 ) 

A 

367 

40300000 


WRITE 16*821 

A 

368 

40400000 


WRITE (6*83) (i*T( l)*xm. YU), 1 = 1,301) 

A 

369 

40500000 


CALL MINMAX IX* 301 *AMIN* AMAX ) 

A 

370 

4060000*0 


CALL MINMAX l Y. 301 * BM IN* 8M AX ) 

A 

371 

40700000 


IF ( AMIN.GT. BMIN) AM IN=BMIN 

A 

372 

40800000 


IF ( AMAX.LT.8MAX) AMAX=8MAX 

A 

373 

40900000 


CALL SCALE BW ( AMIN * AMAX ) 

A 

374 

41000000 


XSAVEINMATH )=AMIN 

A 

375 

41100000 


YSAVEl NMAT+l)=AMIN 

A 

376 

41200000 


SF= I AM AX— AM IN ) / 1 0 • 

A 

377 

41300000 
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XSAVEI NMAT*2)=SF 

A 

378 

41400000 


YSAVE( NMAT *2 )=SF 

A 

379 

41500000 


CALL PLPT ( X SAVE *Y SAVE, NMAT) 

A 

380 

41600000 


CALL INFOPLT ( i*301*X, i*Y* i. AM IN* AMAX* AMIN * AMAX* .5 • 9* XM* 9* YM * 0) 

A 

381 

41700000 


IF IKQ.EQ.4) NKE Y=46 

A 

382 

41800000 

41900000 

* 

COMPUTE RESIDUALS 

A 

383 

420.00000 

42100000 


DO 67 1=2* NMAT 

A 

384 

42200000 


DO 65 J=i * NKR1 

A 

385 

42300000 


J J=NKR 1— J* 1 

A 

386 

42400000 


JI = JJ 

A 

387 

42500000 


IF CTSAVEI l ) .GE.RC JJ ) ) GO TO 66 

A 

388 

42600000 

65 

CONTINUE 

A 

389 

42700000 

66 

IF (JI.GT.NKR) J 1 = NKR 

A 

390 

4280,0000 

67 

Y( l ) = ( ( CS( JI*27*1)*T$AVEC I ) *CS ( J I *27, 2) i *TSAVE( I.) *CS ( J 1*27 * 3 ) )*TSA 

A 

391 

42900000 


1VE { I ) *CSI J 1*27*4) 

A 

392 

43000000 


DO 70 1=2* NMAT 

A 

393 

43100000 


DO 6 8 J= i * NKRS 

A 

394 

43200000 


JJ=NKRS-J*1 

A 

395 

43300000 


JI=JJ 

A 

396 

43400000 


IF ( T S A VE ( 1 1 « GE • 0 ( J J ) ) GO TO 69 

A 

397 

43500000 

68 

CONTI NUE 

A 

398 

43600000 

69 

IF (JI.GT.NKS) JI=NKS 

A 

399 

43700000 


xm»* 1 CSC Jl * 1 ) *TS AVEI I ) *C S ( J I * 2 ) )*TS AVE < I) *CS ( J I, 3 1 )*TSAV E( l)*CS( 

A 

400 

43800000 


1 J I « 4 ) 

A 

401 

43900000 

70 

CONTINUE 

A 

402 

440000(10 


X( 1) = CS( 1,4) 

A 

403 

44100000 


Y( 1)=CS( 23*4) 

A 

404 

44200000 


MAX RES =0- 

A 

405 

44300000 


DO 71 1 = 1 * NMAT 

A 

406 

44400000 


Zin = ABSIXm-XSAVEI I) ) 

A 

407 

4*4500000 


T( I) =ABS(Y I II-YSAVE (I)) 

A 

408 

44600000 


IF (Z( IKGT.MAXRES) MAXR£S=Z(I) 

A 

409 

44700000 


IF < T C I ) • G T« MA XRE S ) MAXRES=TII) 

A 

410 

44800000 

71 

CONTINUE 

A 

41 i 

44900000 


WRITE 16*84) 

A 

412 

45000000 


WRITE (6*8 5) i I • TS AV E C 1 1 * XSAV E( I ) * Y$A VE( I ) * X(I ) , YI l) , Z (I) * T{ I ) * I = I 

A 

413 

45100000 


l .NMAT) 

A 

414 

45200000 


WRITE (6,86) MAX RES 

A 

415 

45300000 


GO TO (39,1,39*1), KO 

A 

416 

45400000 

45500000 

* 


A 

417 

45600000 

45700000 

72 

FORMAT (//9X,1HX, 14X * 1HY * 1 IX* 4HYDQT ) 

A 

418 

45800000 

73 

FORMAT (3E15.6) 

A 

419 

45900000 

74 

FORMAT I IH 1* 10X* , ISH***#******^***^ , / I1X, 15HDAT A FOR X VS T/UX.15 

. A 

420 

46000000 



A 

421 

46100000 

75 

FORMAT ( 40 X * F 7 . 3 * 3 X * F 7 . 3 ) 

A 

422 

46200000 

76 

FORMAT (IX ,1 8HNUMB ER OF CURVES = , 14 * 10X, 18 HNUMBER OF POINTS =*I4/3 

A 

423 

46300000 


12H ENDPOINTS AND JUNCTION POINTS =/ ( 6E20 . 7 ) ) HENDPOI NTSANOJ UNC TIONP 

A 

424 

464000.0 0 


201 NT S= / ( 6E 20 • 7 ) ) 

A 

425 

46500000 

77 

FORMAT (1H0* 12 X, IHT , 19X, 1HY,15X, 10HCDMPUTED Y* 1 OX, 9HRE SI OUAL S ) 

A 

426 

46600000 

78 

FORMAT ( 1H0* 12X, 1HT, 19 X* 1HX*15X* 10HC0MPUTED X , 1 OX, 9HRES I DU ALS ) 

A 

427 

4670,0000 

79 

FORMAT (14*616. 7, 3E20.7) 

A 

428 

46800000 

80 

FORMAT ( /8X* 19HSTANDARD DE VI AT 1 0N=, E 12. 5 ) 

A 

429 

46900000 

81 

FORMAT I IH 1 , 10X, , V5H****^********** * / 11X , 15H0AT A FOR Y VS T/11X,15 

A 

430 

47000000 


IH# a**********#**,//) 

A 

431 

47100000 

82 

FORMAT ( // 3X , IH , 18X • IHT ,20X , 1HX* 20X, 1HY ,/ ) 

A 

432 

47200000 

83 

FORMAT (I5,3E20.8) 

A 

433 

47300000 

84 

FORMAT C.5H1 * 8X , IHT * 9 X* 1HX* 9X,41HY COMPUTED X COMPUTED Y RES X 

A 

434 

47400000 


l RES Y /) 

A 

435 

47500000 

85 

FORMAT (I5*7E10.2) 

A 

436 

47600000 

86 

FORMAT ( // 19 H MAXIMUM RESIDUAL =*E20.7) 

A 

437 

47700000 

87 

FORMAT ( 1H0, 12X, 1HX, 19X* IH Y* 1 5X * 1 QHCQMPUTE D Y * 1 OX, 9HRES I DUALS ) 

A 

438 

47800000 

33 

FORMAT I IH 1 , 10X, , 1 5H***********#***, / 1 IX* 15HDAT A FOR X VS Y/UX,15 

A 

439 

47900000 


IH***#* #**#**##*#,//) 

A 

440 

48000000 


END 

A 

441- 

48100000 
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APPENDIX B - Continued 


Subroutine PLPT. - Subroutine PLPT scales and plots computed values with the 
CalComp POINT routine. The flow diagram for subroutine PLPT is as follows: 



The program listing for subroutine PLPT is as follows: 


SUBROUTINE PLPT CX.YY.NG) 

B 

I 

48200000 

DIMENSION XU). YY 1 1 ) 

B 

2 

48300000 

NP 1=NQ* 1 

B 

3 

48400000 

NP2=NO +2 

B 

4 

48500000 

48600000 

SCALE AND PLOT COMPUTED VALUES WITH POINT ROUTINE 

B 

5 

48700000 

48800000 

XMV=X ( NPl ) 

B 

6 

48900000 

XSF=Xt NP2 I 

B 

7 

49000000 

YMV=YY I NPI ) 

B 

8 

49100000 

VSF= YY I NP2 ) 

B 

9 

49200000 

00 l 1=1* NO 

B 

10 

49300000 

X1=(X( n-XMVJ/XSE 

B 

11 

49400000 

Yl = ( YYCn-YMVl/YSF 

B 

12 

49500000 

CALL POINT (XI. YU 

6 

13 

49600000 

RETURN 

B 

14 

49700000 

END 

B 

15- 

49800000 
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APPENDIX B - Continued 


Subroutine CUSPFIT . - Subroutine CUSPFIT applies the least-squares technique 
described in this report to smooth data using cubic spline functions. The routine obtains 
(1) the value of the cubic spline function and the second derivative at the endpoints and 
junction points, (2) the functional values and values of the first and second derivative for 
a given x, and (3) the coefficients for each segment of the cubic spline function. In the 
flow diagram below, ICODE is a code which specifies the purpose of the current entry into 
the subroutine. ICODE is defined as follows: 

ICODE = 1 Computes second derivatives and functional values at endpoints and 
junction points. 

= 2 Computes functional values for given x. 

= 3 Computes the first derivative for given x. 

= 4 Computes the second derivative for given x. 

= 5 Computes coefficients of cubic equations of the form 
Y = Ax3 + Bx2 + Cx + D 


The flow diagram of subroutine CUSPFIT is as follows: 
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APPENDIX B - Continued 


The program listing for subroutine CUSPFIT is as follows: 




SUBROUTINE CUSPFIT ( ICOOE, C, NKR.L, X , V , R, NMAX.NC , N, A. W.KODE .CS.KLOS 

c 

1 

49900000 



IE) 


c 

2 

50000000 







50100000 

* 


FIELD DEFINITIONS 


C 

3 

50200000 

* 


!CODE=l COMPUTES SECOND DERIVATIVE AND Y VALUES AT R VALUES 


C 

4 

50300000 

* 


IC0QE=2 COMPUTES THE VALUE OF THE FUNCTION FOR GIVEN X VALUES 


c 

5 

50400000 

* 


ICODE=3 COMPUTES THE VALUE OF THE 1ST DERIVATIVE FOR VALUES OF 

X 

c 

6 

50500000 

* 


ICQDE=4 COMPUTES THE VALUE OF THE 2ND DERIVATIVE FOR VALUES OF 

X 

c 

7 

5060.0000 

* 


IC0DE=5 COMPUTES AN EQUATION OF THE FORM Y=AX**3+8X**2+CX+D 


c 

8 

50700000 

* 


XII) AND YU) ARE THE CO-ORDINATES OF THE DATA POINTS 


c 

9 

50800000 

* 


RiJ) IS THE VALUE OF A JUNCTION POINT OR AN END POINT 


c 

10 

50900000 

* 


NMAX IS THE MAXIMUM NUMBER OF POINTS PER CURVE 


c 

11 

51000000 

* 


L0=3*NUMBER OF CURVES+l 


c 

12 

5110000 0 

* 


A IS THE COEFFICIENT MATRIX 


c 

13 

51200000 

* 


C IS THE MATRIX WHICH REPRESENTS THE EQUATION 


c 

14 

51300000 

♦ 


T T T 


c 

15 

51400000 

* 


(A A B MX) (A Y i 


c 

16 

51500000 

* 


( )< ) = I ) 


c 

17 

51600000 

* 


IB 0 UK) I G > 


c 

18 

51700000 

* 


K IS THE L AGRANG IAN MULTIPLIER 


c 

19 

5180000 0 

* 


LIl) ERROR CODE 


c 

20 

51900000 

* 


LU+1) IS THE LENGTH ALONG THE X-AXIS OF THE JTH SEGMENT 


c 

21 

52000000 

* 


NIJ) IS THE NUMBER OF POINTS PER CURVE NINUMBER OF CURVES + 1) = 


C 

22 

52100000 

* 


THE TOTAL NUMBER OF POINTS 


c 

23 

5220000 0 







52300000 



DIMENSION CS (54*41 


C 

24 

524000XX0 



DIMENSION AINMAX.l), CINC.U, Nil). RU), L(l)» XU), YU), IPIVOT 

c 

25 

5250000 0 



un, wui 


c 

26 

52600000 



REAL L 


c 

27 

52700000 



JJ = 0 


c 

28 

5280000 0 



Li n=o. 


c 

29 

52900000 



L0=3*NKR*1 


c 

30 

5300,0000 



IF (KLOSE.EQ.l) LO=LO *2 


c 

31 

5310000 0 



LA-LO 


C 

32 

53200000 



LR=LG*i 


c 

33 

53300000 



LL=L R-NKR 


c 

34 

53400000 



IF (KLOSE.EQ.l) LL=LL-2 


c 

35 

53500000 



NKR!=NKRU 


c 

36 

53600000 



IF MC0DE.EQ.5I GO TO 30 


c 

37 

53700000 



J=2 


c 

38 

53800000 



ISTUNT=0 


c 

39 

53900000 


1 

I COUNT— 0 


c 

40 

54000000 


2 

ICOUNT=IC0UNTM 


c 

41 

54100000 



I 5TUNT= I ST UNT + 1 


c 

42 

54200000 



IF I X! ISTUNTI-RI J) ) 2,3,4 


c 

43 

54300000 


3 

N{ J— i ) = ICOUNT 


c 

44 

54400000 



GO TO 5 


c 

45 

54500000 


4 

Hi J- 1) - 1 COUNT-1 


c 

46 

54600000 



IS rUNT.s I STUNT-1 


c 

47 

54700000 


5 

IF (NIJ-U-NMAX) 7,7,6 


C 

48 

54800000 


6 

JI= J-l 


c 

49 

54900000 



L( 11=2. 


c 

50 

55000000 



PRINT 40, NI JII.JI 


c 

51 

55100000 



GO TO 39 


0 

52 

5520.0000 


7 

J= J*l 


c 

53 

55300000 



IF IJ-NKRU 1,1,8 


C 

54 

55400000 


8 

NINKR*1)=I STUNT 


C 

55 

55500000 



IF I I CODE- 2) 9,LL,24 


c 

56 

55600000 


9 

DO 10 J= 1 • LR 


c 

57 

55700000 



DO 10 1=1, LO 


c 

58 

55800000 



C( I , J 1=0.0 


C 

59 

55900000 


10 

CONTINUE 


c 

60 

56000000 







56100000 
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* 


COMPUTES COEFFICIENT MATRlXU) 

C 

61 

56200000 

56300000 


ll 

IN=Q 

C 

62 

56400000 



IM= 1 

c 

63 

56500000 



n=o 

c 

64 

56600000 



00 20 J=2» NKR1 

c 

65 

56700000 



IK=N< J-l> 

c 

66 

56800000 



IF < IK.EO.O) GO TO 19 

c 

67 

5690000,0 



L(J)=RIJ)-R(H) 

c 

68 

57000000 



IN=IK>IN 

c 

69 

57100000 



DO 12 1= l M * I N 

c 

70 

57200000 



K= I-1M + 1 

c 

71 

57300000 



xi«xt.n-R( Ji 

c 

72 

57400000 



XI Mi = X 1 U-R{ J-1J 

c 

73 

57500000 



A( Kt n={-XI**3/T6.*LI J) )+L(J)*XlMl/6.-L{ Jl**2A6.l 

C 

74 

57600000 



A| K* 2 ) = ( 1 • —X IMl/L ( J ) } 

c 

75 

57700000 



A( K* 31 =(XIM1**3/(6.0*L{J) J-Ll J) *XtMl/6.J 

c 

76 

57800000 



A( K# 41 -XIM 1/LC J i 

c 

77 

57900000 


12 

CONTINUE 

c 

78 

580,00000 



IF (I CODE. EG. 2) GO TO 22 

c 

79 

58100000 

58200000 

* 


WEIGHT POINTS W*A*X=W*Y 

c 

80 

58300000 






58400000 



DO 13 I = I M # IN 

c 

81 

5850.0000 



K= I-IM+l 

c 

82 

58600000 



AIK,lls.A(K.l»*HM.I 

c 

83 

58700000 



At K* 21 = AI K «2 ) *W( I ) 

c 

84 

58800000 



ALK* 31 =A( K*3 ) *W l 1 1 

c 

85 

58900000 



A( K* 41 =A( K » 4 ) *W ( I ) 

c 

86 

59000000 


13 

ym=wcn*Y( n 

c 

07 

59100000 

59200000 

* 


COMPUTES A TRANSPOSE * A 

c 

88 

59300000 

59400000 



KK=0 

c 

89 

59500000 



IL=1 

c 

90 

59600000 



DO 16 K=l,4 

c 

91 

59700000 



KK=K*-l I 

c 

92 

59800000 



DO 15 M= I L ♦ 4 

c 

93 

.59900000 



MM=M+ l I 

C 

94 

60000000 



DO 14 1=1 » I K 

c 

95 

601000,00 



C l KK * M Ml = A II f X 1 * A { I f M I >C I KK , MM I 

c 

96 

60200000 


14 

CONTINUE 

c 

97 

60300000 



C ( MM t KK) =.C IKK.* MM1 

c 

98 

60400000 


15 

CONTINUE 

c 

99 

60500000 



IL= IL* 1 

c 

100 

60600000 


16 

CONTINUE 

c 

101 

60700000 

60800000 

* 


COMPUTES A TRANSPOSE * Y 

c 

102 

60900000 

61000000 



DO 17 M=l ♦ 4 

c 

103 

61100000 



MM=M+I I 

c 

104 

61200000 



DO 17 1=1 M * IN 

c 

105 

61300000 



K=I-CM + 1 

c 

106 

61400000 



C{MM,LR1=A<K,MI*Y ( I)+CIMM,LRI 

c 

107 

61500000 


17 

CONTINUE 

c 

108 

61600000 



IF i J.E0.2) GO TO 19 

c 

109 

61700000 

61800000 

* 


COMPUTES CONDITIONS FORCING THE FIRST DERIVATIVE TO BE CONTINUOUS 

c 

110 

61900000 

62000000 



J=J-i 

c 

111 

62100000 



LL J=LL+ J-1 

c 

112 

62200000 



K= I 1-2 

c 

113 

62300000 



C ( LL J* l + K) =L(J)/6» 

c 

114 

62400000 



CC LL J* 2+K J=— 1 */L C J ) 

c 

115 

62500000 



ctLLJ*3frK)=(L( ji autin/3. 

c 

116 

62600000 



c(LLJr4^Ki-(Li j)>L{jtin/(Uj>n*L( j.i > 

c 

117 

62700000 



Cl LL J ♦ 5+K ) =L I J #■ 1 ) / 6« 

c 

118 

62800000 



CILLJ.6+K)=-1./LU + 1> 

c 

119 

62900000 



C< 1+K* LL J ) =C ( LL J * li-K) 

c 

120 

63000000 



C l 2 + K* LL J I =C (LL J* 2+K) 

c 

121 

63100000 
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* 


* 


18 

19 

20 
21 


* 


22 


23 

24 


* 


25 


Ct 3 FK« LL J I =C (LLJ.3FK1 

c 

122 

63200000 

CI4fK.LLJ>=C(LLJ,4fK) 

c 

123 

63300000 

C( 5fK« LL J 1 =C( LL J , 5fk) 

c 

124 

63400000 

C16FK.LLJI =C (LL J, 6FK ) 

c 

125 

63500000 

IF (KLOSE.EO.OI GO TO 18 

c 

126 

63600000 

IF I J.NE.NKRl-l) GO TO 18 

c 

127 

63700000 

63800000 

FORCE EQUALITY OF THE FIRST DERIVATIVES AT THE ENDPOINTS 

c 

128 

63900000 

64000000 

LLJ=LL JfI 

c 

129 

64100000 

CILLJ, 3fK)=-L(Jf1)/6. 

c 

130 

64200000 

CILLJ. 4FKl=L./LtjFl) 

c 

131 

64300000 

CILLJ. 5fK)=-LLJF1)/3. 

c 

132 

64400000 

C( UL J* 6*K 1 =— l • /L I J F 1 1 

c 

133 

64500000 

CILLJ. l»=-LI2)/3. 

c 

134 

64600000 

C(LLJ»2)=-1./L(2) 

c 

135 

64700000 

CILLJ, 3)=-L(2J/6. 

c 

136 

64800000 

CILLJ, 41=1. /LI2I 

c 

137 

64900000 

Cl 3FK.LL J ) =C I LL J , 3 FK I 

c 

138 

65000000 

C I 4 F K . LL J I =C 1 LL J , 4 FK 1 

c 

139 

65100000 

Cl 5FK, LL J )=C (LL J , 5 FK) 

c 

140 

65200000 

Cl 6FK,LLJJ=CILLJ,6FK) 

c 

141 

65300000 

Cl 1, LL J ) = C ( LLJ, 1 ) 

c 

142 

65400000 

C(2»LLJ)=C(LLJ«2) 

c 

143 

65500000 

Cl 3.LLJ1=CILLJ.3) 

c 

144 

65600000 

C(4,LLJ)=CILLJ,4I 

c 

145 

6570.0000 

65800000 

FORCE EQUALITY OF THE FUNCTIONAL VALUES AT THE ENDPOINTS 

c 

146 

65900000 

66000000 

LL J=LL JFl 

c 

147 

66100000 

CILLJ, 21=1. 

, c 

148 

66200000 

Cl 2, LLJ 1 = 1 , 

c 

149 

66300000 

CILLJ, 6FK)=-1. 

c 

150 

66400000 

Cl 6FK, LLJ) =-l • 

c 

151 

66500000 

J= J F l 

c 

152 

66600000 

1 1 = I IF2 

c 

153 

66700000 

IM= I NF 1 

c 

154 

66800000 

CONTINUE 

c 

155 

66900000 

IF IICODE.EQ.2) GO TO 39 

c 

156 

67000000 

CALL SIMEO I C,L A ,C 1 1 ,LR) , 1 , DET , L, NC » I SCALE ) 

c 

157 

67100000 

IF (ICODE.EQ.5J GO TO 32 

c 

158 

6720.0000 

GO TO 39 

c 

159 

67300000 

67400000 

COMPUTES VALUES OF THE FUNCTION 

c 

160 

67500000 

67600000 

00 23 I = IM , I N 

c 

161 

67700000 

K= I-IMFl 

c 

162 

67800000 

YI n=C( lFl I,LR ) *A(K, 1IFC I2FII ,LRI*A1K,2) fCI3fII , LR) *A( K, 3) FC 1 4.F 1 1 , 

c 

163 

67900000 

1LR ) FAIR, 4) 

c 

164 

68000000 

CONTINUE 

c 

165 

68100000 

GO TO 19 

c 

166 

68200000 

IN = 0 

G 

167 

68300000 

IM=1 

c 

168 

68400000 

IF ( IC0DE.EQ.4) GO TO 27 

c 

169 

68500000 

68600000 

COMPUTES VALUES OF THE FIRST DERIVATIVE 

c 

170 

68700000 

68800000 

DO 26 J=2 » NKR1 

c 

171 

68900000 

LI J)=RUI-R( J-1J 

c 

172 

69000000 

IN=NI J-1IFIN 

c 

173 

69100000 

SLJ=C( 2*J-l,LRI 

c 

174 

69200000 

SL J1=C 1 2* J-3 ,LR) 

c 

175 

69300000 

YLJ=C(24J.LR) 

c 

176 

69400000 

YL J1=C l 2* J— 2 »LR) 

c 

177 

69500000 

DO 25 I=IM, IN 

c 

178 

69600000 

Y( 1 1 =SL Jl* l- { Xl I )-R€ J) 1 **2/1 2. *L< J) I FL I J 1/6. I+SLJ*IIXI IL-R(J-I) ) ** 

c 

179 

69700000 

12/ (2. *LIJI I-LIJ 1 / 6.1 FJ YLJ-YLJn/LUl 

c 

180 

69808000 

CONTINUE 

c 

181 

6990000.0 

IM= I NF 1 

c 

182 

70000000 
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26 

CONTINUE 

C 

183 

70100000 



GO ro 39 

C 

184 

70200000 






703000.00 

* 


COMPUTES VALUES OF THE SECOND DERIVATIVE 

c 

185 

70400000 






70500000 


27 

on 29 J=2 , NKR 1 

c 

186 

70600000 



SL J=C ( 2* J- L , LR) 

c 

187 

70700000 



SLJl=C(2*J-3»LRJ 

c 

188 

70800000 



YL J=C ( 2* J » LR J 

c 

189 

70900000 



YLJL=CI2*J-2,LR) 

c 

190 

71000000 



1N=N( J-1M-IN 

c 

191 

71100000 



LI J) = R( JJ-RI J-1.1 

c 

192 

7120.0000 



DO 28 1= I M * I N 

c 

193 

71300000 



YM)=SLJ*I (XCl)-RI J-in/U JII +SLJ 1* < ( RI J l-XI Ilj/UJl) 

c 

194 

71400000 


28 

CONTINUE 

c 

195 

71500000 



I M= l N+ 1 

c 

196 

71600000 


29 

CONTINUE 

c 

197 

71700000 



GO TO 39 

c 

198 

71800000 






71900000 

* 


COMPUTES THE EQUATIONS OF THE FORM AX **3+BX**2+CX*D 

c 

199 

72000000 






72100000 


30 

LR1=LR 

c 

200 

72200000 



DO 37 J=2 , NKRl 

c 

201 

7230.0000 



DO 31 IL=1 ,2 

c 

202 

72400000 



C( 2*IL-1, 2 1=2. 

c 

203 

72500000 



C{ 2* I L— 1 ,31=0. 

c 

204 

72600000 



C( 2'* l L- l , 4 1 = 0* 

c 

205 

72700000 



Cl 2* IL. 41 = 1. 

c 

206 

72800000 



CC 2*IL-i,i J = 6.*R( J+IL-2J 

c 

207 

72900000 



DD 31 1=1,3 

c 

208 

73000000 



C( 1*51 =C t 2* J-4+ I » LR ) 

C 

209 

73100000 



C( 2* IL * LI = I R ( U* I L-2 ) 1 **f 4— I ) 

c 

210 

73200000 


31 

CONTINUE 

c 

211 

73300000 



C(4,5)=C(2*J,LR> 

c 

212 

73400000 



LA=4 

c 

213 

73500000 



LR=5 

c 

214 

73600000 



GO TO 21 

c 

215 

73700000 


32 

LR=LR1 

c 

216 

73800000 



J 1=4-1 

c 

217 

73900000 



PRINT 41, JI,(C( 1,5) •1 = 1,4) 

c 

218 

74000000 






741000.00 

* 


SAVE COEFFICIENTS 

c 

219 

74200000 






74300000 



IF (K0DE.E0.4J GO TO 3o 

c 

220 

74400000 



IF ( KODE.NE. 0) GO TO 34 

c 

221 

74500000 



DO 33 1=1,4 

c 

222 

74600000 


33 

CS ( J I • I }=C ( I ,5) 

c 

223 

74700000 



GO TO 36 

c 

224 

74800000 


34 

DO 35 1=1,4 

c 

225 

74900000 


35 

CS< JI*27.I)=C<I,5) 

c 

226 

75000000 


36 

CONTINUE 

c 

227 

75100000 



NR=LO— NKRl +J 

c 

228 

75200000 



DO 37 1=1,4 

c 

229 

75300000 



Cl NR, 1 1=CI 1,5) 

c 

230 

75400000 


37 

CONTINUE 

c 

231 

75500000 



DO 38 J=l , NKR 

c 

232 

75600000 



NR=LR-NKR1 +J 

c 

233 

75700000 



DO 38 1=1,4 

c 

234 

75800000 



Cl J, I)=C(NR, II 

c 

235 

75900,000 


38 

CONTINUE 

c 

236 

76000000 


39 

RETURN 

c 

237 

76100000 






76200000 

* 



c 

238 

76300000 






76400000 


40 

FORMAT ( 36H ERROR COND IT ION — NMAX I S TOO SMALL. /UH THERE ARE ,15, 

c 

239 

76500000 



1 17H POINTS ON CURVE ,12, IH.) . 

c 

240 

76600000 


41 

FORMAT </7H CURVE , I2/4X.F23.9, 7H< X**3M-,E 23*9, 7H{ X**2 U ,F23. 9,4H( 

c 

241 

76700000 



iXHv,F23.9l 

c 

242 

76800000 



END 

c 

243- 

76900000 
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Subroutine SUP. - Subroutine SUP, used in cases selecting the parametric version, 
computes the parametric variable t and sets up the arrays to be used in subroutine 
CUSPFIT. The flow diagram for subroutine SUP is as follows: 


SUBROUTINE 

SUP 


/' 




Compute 

parametric 

variable 


r 

i 


▼ 

! 

Move arrays to 
make t independent 
variable and x 
I dependent 

I variable 



The program listing for subroutine SUP is as follows: 



SUBROUTINE SUP { X * Y, T , Z . XMAXT , NMAT) 

D 

i 

77000000 

77100000 


COMPUTE PARAMETRIC VARIABLE AND SET UP ARRAYS 

D 

2 

77200000 

77300000 


DIMENSION XI 31 0 ) » YI310I, T( 3 i Qi • ll 310) 

D 

3 

77400000 


T( 1 ) = X MAXT =0 * 

D 

4 

77500000 


00 l 1=2, NMAT 

0 

5 

77600000 

1 

t 1 1 )= ( x 1 1 j — x ( i — i i ) ** 2 * i y i i )-yi i-in**2+TM-ii 

D 

6 

77700000 


IF (Tin. GT. XMAXT) XMAXT = TU) 

0 

7 

77800000 


DO 2 1=1, NMAT 

D 

8 

77900000 


71 I)=Y( u 

D 

9 

78000000 


YC I 1 — XX II 

D 

10 

78100000 

2 

xi n=rm 

0 

11 

78200000 


RETURN 

D 

12 

78300000 


END 

D 

13- 

78400000 
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Subroutine MINMAX . - Subroutine MINMAX finds the minimum and maximum values 
of data to be plotted when using the parametric option of the program. The flow diagram 
for subroutine MINMAX is as follows: 




The program listing for subroutine MINMAX is as follows: 


SUBROUTINE MINMAX I A. N. AMI N* AM AX) 

* TO FIND MINIMUM AND MAXIMUM VALUES IN 

* C A= AN ARRAY . USER MUST DIMENSION A TO FIT HIS ARRAY 

* C N- NO OF VALUES IN THE ARRAY 

DIMENSION At 310) 

AM l N= I • 0E2 0 
AM AX= I • 0E- 20 
DO I 1 = 1. N 

IF I At I ) . L T* AMIN I AMIN=AU) 

IF t At D.GT.AMAX) AM AX- A ( I ) 

1 CONTINUE 
RETURN 
END 


E 1 78500000 

78600000 
E 2 78700000 

E 3 78800000 

E * 78900000 

79000000 
E 5 79100000 

E 6 79200000 

E 7 79300000 

E 8 79*00000 

E 9 79500000 

E 10 79600000 

E 11 79700000 

E 12 79800000 

E 13- 79900000 
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Subroutine SCALEBW . - Subroutine SCALEBW scales the values for plotting in cases 
using the parametric option. The flow diagram for subroutine SCALEBW is as follows: 



The program listing for subroutine SCALEBW is as follows: 


SUBROUTINE SCALEBW I YMIN# YM AX ) 

F 

1 

80000000 

80100000 

SCALE FOR PLOTTING 

F 

2 

8020,0000 

80300000 

DIMENSION FAC <31 

F 

3 

80400000 

FACT! 1=1 .0 

F 

4 

80500000 

FACI 21-2.0 

F 

5 

80600000 

FAC I 3 1=5.0 

F 

6 

80700000 

FAK = 1 • OE-8 

F 

7 

80800000 

YM INS= YMIN 

F 

8 

80900000 

YMAXS= YMAX 

F 

9 

81000000 

A= ( YMAX— YM IN ) / 10 .0 

F 

10 

81100000 

CONTINUE 

F 

11 

81200000 

DO 3 J=l* 16 

F 

12 

81300000 

FAK=10.0*FAK 

F 

13 

81400000 

DO 2 1=1# 3 

F 

14 

81500000 

C=FAK*FAC( I) 

F 

15 

81600000 

IF (A.GT.C) GO TO 2 

F 

16 

81700000 

B=AM0D{ YMIN.C) 

F 

17 

81800000 

IF ( YMIN.GT.0#0J YMI N=YM IN— B 

F 

18 

81900000 

IF ( YMIN. LT. CKO) Y M I N= Y M I N— I C + 8 ) 

F 

19 

82000000 

YMAX=YMIN+10.0*C 

F 

20 

82100000 

IF JYMINS.LT. YMIN) GO TO 4 

F 

21 

82200000 

IF ( YMAXS.GT.YMAXi GO TO 4 

F 

22 

82300000 

RETURN 

F 

23 

82400000 

CONTINUE 

F 

24 

82500000 

CONTINUE 

F 

25 

82600000 

PRINT 5, A # YMI N# YMAX*C 

F 

26 

82700000 

RETURN 

F 

27 

82800000 

A={ YMAXS-YMINSI/5.0 

F 

28 

82900000 

GO TO l 

F 

29 

83000000 

83100000 


F 

30 

83200000 

83300000 

FORMAT ( /40H NEED TO ALTER SUBROUTINE SCALEBW A=E15. 8# 2X# 5HYMI 

F 

31 

83400000 

1N=E15 • 8#2X #5HYMAX=E15.8# 2X »5M$$$$$/3X#2HC=E15*8/) 

F 

32 

83500000 

END 

F 

33- 

83600000 
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Langley Library Subroutine SIMEQ 

Language: FORTRAN 

Purpose : SIMEQ solves the matrix equation AX = B where A is a square coefficient matrix and B is a 
matrix of constant vectors. The solution to a set of simultaneous equations and the determinant may be 
obtained. If the user wants the determinant only, use DETEV for savings in time and storage. 

Use: CALL SIMEQ (A, N, B, M, DETERM, IPIVOT, NMAX, ISCALE) 

A A two-dimensional array of the coefficients. 

N The order of A; 1 .£ N £ NMAX. 

B A two-dimensional array of the constant vectors B. On return to calling program, X is 

stored in B. 

M The number of column vectors in B. 

DETERM Gives the value of the determinant by the following formula: 

DET(A) = 10 100 ISCALE (DETERM) 

IPIVOT A one -dimensional array of temporary storage used by the routine. 

NMAX The maximum order of A as stated in dimension statement of calling program. 

ISCALE A scale factor computed by subroutine to keep results of computation within the floating- 

point word size of the computer. 

Restrictions: Arrays A, B, and IPIVOT are dimensioned with variable dimensions in the subroutine. The 
maximum size of these arrays must be specified in a DIMENSION statement of the calling program as: 

A (NMAX, NMAX), B (NMAX, M), IPIVOT (NMAX). The original matrices, A and B, are destroyed. 

They must be saved by the user if there is further need for them. The determinant is set to zero for 
a singular matrix. 

Method: Jordan's method is used through a succession of elementary transformations: l n , l n -i? . . l \ * 

If these transformations are applied to a matrix B of constant vectors, the result is X where AX = B. 
Each transformation is selected so that the largest element is used in the pivotal position. 

Accuracy: Total pivotal strategy is used to minimize the rounding errors; however, the accuracy of the 
final results depends upon how well -conditioned the original matrix is. 

Reference: (a) Fox, L.: An Introduction to Numerical Linear Algebra. Oxford Univ. Press, c. 1965. 

Storage: 4320 locations. 

Subroutine date : August 1, 1968. 
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Langley Library Subroutine FTLUP 


Language : FORTRAN 

Purpose : Computes y = F(x) from a table of values using first- or second-order interpolation. 

An option to give y a constant value for any x is also provided. 

Use : CALL FTLUP(X, Y, M, N, VARI, VAKD) 

X The name of the independent variable x„ 

Y The name of the dependent variable y = F(x). 

M The order of interpolation (an integer) 

M = 0 for y a constant. VARD® corresponds to VARI(I) for 
1=1,2,. . ., N. For M = 0 or N i 1, y = F(VARI(1)) for any value of x. 

The program extrapolates. 

Ms 1 or 2. First or second order if VARI is strictly Increasing (not equal). 

M = -1 or -2. First or second order if VARI is strictly decreasing (not equal). 

N The number of points in the table (an integer). 

VARI The name of a one -dimensional array which contains the N values of the independent variable. 

VARD The name of a one -dimensional array which contains the N values of the dependent variable. 

Restrictions : All the numbers must be floating point. The values of the independent variable x in the 
table must be strictly increasing or strictly decreasing. The following arrays must be dimensioned by 
the calling program as indicated; VARI(N), VARD(N) . 

Accuracy : A function of the order of interpolation used. 

References : (a) Nielsen, Kaj L.: Methods in Numerical Analysis. The Macmillan Co., c.1956, pp. 87-91. 

(b) Milne, William Edmund: Numerical Calculus. Princeton Univ. Press, c.1949, pp. 69-73. 

Storage : 430g locations. 

Error condition : If the VARI values are not in order, the subroutine will print TABLE BELOW OUT OF 
ORDER FOR FTLUP AT POSITION xxx TABLE IS STORED IN LOCATION xxxxxx (absolute). It then 
prints the contents of VARI and VARD, and STOPS the program. 

Subroutine date : September 12, 1969. 
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Usage 

The program D3670 is run on the Control Data series 6000 computer under the 
scope 3.0 operating system. The storage required for a batch run is approximately 
70000g locations and for an on-line CRT rim is approximately 75000g locations. Cases 
implementing the on-line CRT capability use the CDC 250 CRT console. Instructions are 
written into the program and displayed on the screen to inform the on-line CRT user of 
options available to him at various points in the program. The user conveys his selec- 
tions through the use of the function keyboard and the typewriter keyboard at the CRT con- 
sole (see fig. 3). If the curve fit obtained is not satisfactory, an on-line CRT user can 
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APPENDIX B - Continued 



and possibly the number of curves (m 


1) 


in the spline approximation and go back to recompute the spline function. When running 
a parametric CRT case, the user can go back at any time to recompute S^ x (t) or S^ y (t) 
so that these approximations might be improved. This can be desirable after displaying 
the curve on an x-y coordinate system with points computed as t is incremented through 
its range. 


When a case is run as a batch job, the options mentioned above are not available so 
that the job is completed with the values originally input. 


Input Description 

Input is standard CDC FORTRAN NAMELIST. There is an option available allowing 
the coordinates of the input points to be read in a format so that cards punched by another 
program could be directly input into this program. If a batch parametric case is being 
run, the first input values for the joints are taken as the joints for the spline function 
S a (t). A second set of input is read so that the joints may be changed for computing 
SAy(t). If a CRT parametric case is being run, the second read could be bypassed if 
the same joints used in the computation for SA„(t) are desired for the computation of 
SAy(t). This is done according to instructions from the CRT console. 

To simplify the necessary input in a standard nonparametric case, the value of the 
abscissa of the first point to be smoothed is subtracted from all other abscissas before 
the spline functions are computed. This enables the user to input joints beginning with 
zero and stepping up to the value of the abscissa of the last point minus the value of the 
abscissa of the first point. 

The NAMELIST input data, listed under $NAM1 are given as follows: 

X array of abscissas of points to be smoothed 

YY array of ordinates of points to be smoothed 

NKR number of cubics to be used to fit the spline function 

(DEFAULT = 1) 

R array of (NKR + 1) values for the endpoints and joints between curves in 

the spline function 

(DEFAULT: R(l) = 0.0, R(NKR + 1) = X(NMAT) - X(l)) 

NMAT number of points to be smoothed 
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IFLAG 

= 0 if coordinates of points are input with a format 
= 1 if coordinates are input with NAMELIST 
(DEFAULT - 1) 

W 

array of NMAT values of weights for input points 
(DEFAULT = 1.0, i = 1, 2, . . ., NMAT) 

K0 

= 1 for standard CRT version 
= 2 for standard batch version 
= 3 for parametric CRT version 
= 4 for parametric batch version 
(DEFAULT = 1) 

KL0SE 

= 0 for arbitrary curve 
= 1 for closed curve parametric case 
(DEFAULT = 0) 


Output Description 

Output is in the form of plotted curves and printed data. For a standard case, the 
input data points and the computed spline function are plotted on an x-y grid. Small ver- 
tical bars are plotted to indicate the junction points in the spline function. For an on-line 
case, these joints may then be manipulated to obtain a satisfactory curve fit. The printed 
data for a standard case consist of a listing of the number of points, the number of curves, 
the endpoints and junction points, and a table giving the input y values, the computed 
spline and the residuals at the input x values. This is followed by the residual standard 
deviation and a listing of the coefficients of the computed cubic curves making up the 
spline function. 

For a parametric case, similar data is plotted on a t-x grid and printed with this 
input and computed x as a function of t. This is followed by a plot on a t-y grid and 
printed output of y as a function of t. The input points are then plotted on an x-y grid 
with the curve determined by incrementing through values of t and plotting S^y(SA x (t)). 
A listing of the value of t for the input x and y is also printed with x, y, S^ x (t), 
SAy( s A x (t)), and the residuals. 
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APPENDIX C 


EXAMPLE APPLICATIONS OF ALGORITHMS AND COMPUTER PROGRAM 


This appendix describes three example applications of the use of the algorithms and 
corresponding computer program. These examples demonstrate the explicit techniques 
and the parametric technique with and without the closed-curve option. The data sets of 
the applications are chosen for their demonstrative character rather than their relative - 
ness to particular engineering problems. The input data are shown as they would be 
written for submittal to the computer, and the output data are shown as they would appear 
on the CRT and output listing. 


Example Applications 

Case 1 Explicit Algorithm 

Case 1 is a standard case, the input is as follows: 


$NAM1 

X = 0.0, 0.1E+00, 0.2E+00, Q.3E+00, 0.4E+0C, 0.5E+00, O.fcE+OO, 

0.7E+00, 0. 8E + 0C » 0.9E+0C, 0.1E+01, 0.11E+01, 0.12E+01* 

0. 13001, 0.14E+01, 0.15E+01, 0.16E+01, 0.17E+01, 0.18E+01, 

0 • 19E+01 » 0.2E+01, 0.21E+C1, 0.22E+01, 0.23E+01, 0.24E+01, 

0 . 25E+01 , 0.26E+01* 0.27E+01* 0.28E+01, 0.29E+01, 0.3E+01, 

Q.31E+01* 0.32E+01, 0.33E+01, 0.34E+01, 0.35E+01, 0.36E+01, 

0.37E+01, 0.38E+01, 0.39E+01, 0.4E + 01, 0.41E+01, 0. 42E+01, 

0.43E+01* Q.44E+01, 0.45E+C1, 0.46E+01, 0.47E+01, 0.486+01, 

0.49E+01, 0.5E+01, 0. 516+01, 0.52E+01, 0.53E+Q1, 0. 546+01, 

0 • 55E + 01 » 0.56E+01, 0.57E+C1, 0.58E+01, 0.59E+01, 0.6E+01, 

0.61E+01* 0.62E+01, 0.63E+C1, 0.64E+01, 0.65E+01, 0.66E+01, 

0. 67E+01 , 0.68E + 01, 0.69E+01, 0.7E+01, 0.71E+01, 0.72E+01, 

0.73E+01, 0.74E+01, 0.75E+Q1, 0.76E+C1, 0. 77E+01* 0.78E+01, 

0. 79E+01 , 0.8E + 01, 0.81E+01, 0.826+01, 0.83E+01* 0.84E+01, 

0.85E+01, O.86E+01, 0.87E+C1, 0.88E+01, 0.89E+01, 0.9E+01, 

0. 916+01, 0.92E+01, 0.93E+01, 0.94E+01, 

YY = — 0.311E+01 , -0.32E+01, -0.309E+01, -0.315E+01, -0.306E+01, 

-0.285E+01, -0.268E+01, -0.256E+01, -0.232E+01, -0.244E+01, 
-0.252E+01, -0. 239E+01 , -C.24E+01, -0.256E+01* -0.238E+01, 

-0 • 225E +01 , -0.229E+01, -0.225E+01* -0.209E+01, -0.216E+01, 
-0.188E+01, -0.213E+01, -0. 216+01, -0.161E+01, -0.165E+01, 

— 0. 191994E+01, -0. 166E + 01 , -0.192E+01* -0.244E+01* -0.211E+01, 

-0 .2216+0 1 , -0.222E+01, -0.245E+01, -Q.245E+01, -0.26E+01, 

-0. 263E+01 » -0. 281E+01 , -0.323E+01, -0.347E+01, -0.345E+01, 

— 0 .364E+Q1 , -0.361E+01, -0.35E+01, -0.316E+01, -0.315E+01 , 

-0.311E+01, -0.3056 + 01, -0.265E+01, -0.301E+01, -0.317E+01, 

-0.32E+01, -0.323E+01, -0.353E+01, -0.395E+01, -0.376E+01, 

— 0. 379E+01 , -0.384E+01, -0.399E+01, -0.378E+01, -0.378E+01, 
-0.328E+01, -0 . 352E+01 , -0. 3446+01, -0.305E+01, -0. 302E+01, 

-0. 261E+01 , -0.287E+01, -0.262E+01, -0.219E+01, -0.288E+01, 

-0.34E+01, -0.3466+01, -0.381E+01, -0.408E+01, -0.436+01, -0.434E+01, 
— 0.416E+01 , -0. 502E+01 » -0.412E+01, -0.337E+01, -0.291E+0I, 
-0.2C8E+01, -0 .1296+01, -0.41E+0C* 0.11E+00, 0.9E-01, 0.63E+00, 

O.IE+Oi, 0.43E+00, 0.26E+00, 0.91E+00, 0.167E+01, 0.328E+01, 

0.5C3E+01, 0.715E+01, 
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NKR - 2, 

R = 0.0. G.35E+01, 0.94E+0 1 , 


NMAT 

.as 

95, 



I FLAG 


If 



W 

s 

0 .IE+01 , 

O.IE+OI, 

0. 



0. IE+01, 

O.IE+OI, 

0. 



O.IE+OI, 

O.IE+OI, 

0. 



o.ie+oi. 

0. IE+01, 

0. 



0. IE+01, 

O.IE+OI, 

0. 



0. IE+01, 

O.IE+OI, 

0. 



0 • IE + 01 , 

O.IE+OI, 

0. 



0. IE+01, 

O.IE+OI, 

0. 



O.IE+OI, 

O.IE+OI, 

0. 



0. IE+01, 

O.IE+OI, 

0. 



0. IE+01, 

0. IE + 01 , 

0. 



0. IE+01, 

O.IE+OI, 

0. 



0 • IE >.01:9 

0. IE+01, 

0. 



O.IE+OI , 

O.IE+OI. 

0. 

KC 

= 

2, 



KLOSE 

- 

Ot 



$ END 






IE+01 , 

0. IE + 01, 

0. IE+01, 

0. IE+01 , 

IE + 01. 

0. IE+01, 

0. IE+01, 

0. IE+01, 

IE + 01 , 

0. IE+01, 

0. IE+01, 

0. IE+01, 

IE +01 , 

0. IE+01, 

0. IE+01, 

0. IE+01, 

IE+01, 

0. IE + 01, 

0. IE+01, 

0. IE+01, 

IE+01 , 

0. IE+01, 

0. IE+01, 

0. IE+01, 

IE+01, 

0. IE+01, 

0. IE+01, 

0. IE+01, 

IE+01, 

0. IE+01, 

0. IE+01, 

0. IE+01, 

IE+01, 

0. IE+01, 

0. IE+01, 

O.IE+OI, 

IE+01, 

0. IE+01, 

0. IE+01, 

0. IE+01, 

IE+01, 

O.IE+OI, 

0. IE+01, 

0. IE+01, 

IE+01, 

0. IE+01, 

0. IE+01, 

0. IE+01, 

IE+01, 

0. IE+01, 

0. IE+01, 

0. IE+01, 

IE+01, 

0. IE+01, 




In this case, it is desirable to find a curve fit which will smooth the 
The output is as follows: 


0. IE+01, 
0. IE+01, 
0. IE+01, 
0. IE+01, 
O.IE+OI, 
0. IE+01, 
0. IE+01 , 
0* IE + 01 y 

o.ie+oi, 

0. IE+01, 
0. IE+01 , 
0. IE+01, 
0. IE+01, 


input data. 
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NUMBER OF CURVES = 2 NUMBER OF POINTS ? 95 

ENDPOINTS AND JUNCTION POINTS = 

0. 3. 5000000E+00 9.4000 OOOE+OO 


APPENDIX C - Continued 
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The output from this improved fit is as follows: 





NUMBER Of CURVES * 16 NUMBER Of POINTS * 95 

ENDPOINTS AND JUNCTIUN POINTS * 

0. 1 .40000006+00 2.00000006+00 2.50000006+00 3.00000006+00 4.00000006+00 

4.5000000t+00 5.00000006+00 5. 50000006+00 6.00000006+00 6. 5000000E+00 7. OOOOOOOE+OO 
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47941 8546 ( X**3 ) *■ 172 ,049992849< X**2)t - 1758. 3709519121 X)> 5884.647612931 



APPENDIX C - Continued 


Case 2 


Case 2 is a parametric algorithm, closed curve. The input with junction values 
for SA x (t) is as follows: 

SNAMl 

X = -0. 131436+01, -0.140676+01, -0.14576*01, -0.157256*01, -0.168976*01, 

-0.17946*01, -0.191436*01, -0. 20442E *0 1 , -0.215216*01, -0.223766*01 , 

-0. 24346+01, -0.255766+01 , -0.26866+01, -0. 283 53E+01 , -0.29571E+01, 

-0.333216+01, -0.316226+01, -0.324396+01, -0.33O1E+O1, -0.333616*01 , 

-0.333746+01, -0. 332296+01, -0.32796+01, -0.323746+01, -0. 318 526+0 1 , 

-0.313246+01, -0. 305426+01, -0.297716+01, -0.2926E+01, -0. 2 86296+01 , 

-0.276636+01 , -0. 266176+01, -0.256986+01, -0.24946^+01, -O.24114E+0 1, 

-0. 233486+01, -0.225346+01, -0.216796+01, -0. 20568E +0 1, -0.197956+01, 

-0.188036+01, -0. 130246+01, -0.174666+01, -0.167966 +01, -0. 163 1 1 E+0 1 , 

-0.159726+01, -0.156676+01, -0.151566+01, -0.14786+01, -0.143576+01 , 

- 0. 1 36 076 + 01 , -0. 136556*01, -0.131246 + 01, -0. 123516+01, -0.115786+01, 

-0. 103046+01, -0.100316+01, -C. 92526+00, -0.92196+00, -0.69986+00, 

-0. 5958E + 00, —0. 4733E + 0Q , -0.34676 + 00, -0. 25146 + 00, -0. 14176 + 00, 

-0.2916-01, 0.9426-01, 0.23366+00, 0.35946+00, 0.48496+00, 

0.60456+00, C. 73256+00, 0. 95496+00, 0.96926+00, 0.106836+01, 

O.113O3E+01, 0.120536+01, 0.128046+01, 0. 135546+01, 0.143046+01, 

C. 1489E + Q1 , 0.155026+01 , 0. 159996 + 01, 0.16194E+C1, 0.166 846+01, 

0.171366*01, 0. 175296 + 01, 0.181056 + C1, 0. 136316 + 01, 0. 1929 IE +01, 

0.197866+01, 0. 2O503E+O1, 0.212186+01, 0.220456+01, 0.22802E+0 1, 

0. 240056+01, 0.248366+01, 0.2591 36+01, 0.276+01 , 0.278626+01, 

0.287476+01, G. 30052 E +01, 0. 311236+01, 0.319656+01, 0.326216+01, 

C.33141E+01, 0.33142E+01, 0.331166+01, 0.32931E+01, 0.326716+01, 

0. 32446E + 01 , 0. 31976E + 01 , G.31609E + 01, Q.30936E+01, 0. 300686+01, 

0.290536+01, 0. 278796+01, 0. 267266+01, 0.254566*01, 0.242766+01, 

0.230746+01, 0. 222056+01, 0.211776+C1, 0. 201426+01, 0.191316+01, 

0*18066+01, 0.17079R+C1 , 0. 16002 E +01 , 0.149916+01, 0.142566+01, 

0.137366+01, 0. 133776+01, 0. 126536+01, 0. 119796+01, 0.111496+01, 

0. 1C3E + 0 i , 0. 94 84E + 00, 0.8 3266 + 00, 0.7162 6+00, 0.58636+00, 

0.43256+00, 0. 31266+00, 0.2216+00, 0. 1028E+OO, -0.276E-O1, 

-0.14046 + 00, -0.27736+00, - 0. 40046 + 00, -0.51746 + 00, -0.63026 + 00, 

-G.74°9fc+00, -u. 84166+00, -0.93176+00, -0. 102 G7E+0 1 , -0.111466+01, 

-C. 122546*01, -0.13143 E +01, 

YY = -0.129836+01, -0. 12234E+01 , -0. 1 1862E+0 1 , -0. 11143E +01 » -0. 10539E+0 1, 

-0.100276 + 01, — 0. 94 336 +00 , -0.88266 + 00, -0.83626 + 00, -0.78246+00, 

-0. 72786+00, -0.68486+00, -0.63966+00, -0.58396+00, -0. 53546+00, 

- 0* 4818b +0 0, -0.4443E+OO, -0. 39286 + 00, -0. 32936+00, -0.22886 + 00, 

-G. 15036+00, -0.1273E+OO, -0.6516-01, -0.2576-01, 0.1496-01, 

0. 5186-01 , 0.10536+00, 0. 15446+00, 0. 18596+00, 0.2219E+00, 

C.2733E+00, O.3265E+00, 0.373E+00, C. 41 146+00, G. 45346+00, 

0.49^76 + 00, 0.53646 +00, 0.57926 + 00, 0.63646+00, 0.67786 +00, 

O.7327E+0O, 0.77646+00, 0.80956+0 0, 0.84586+00, 0.87156+00, 

0.89666*0 0, 0.9466E + 00 , 0.10226E + 01, 0. 1O791E+01 , O.11339E+01, 

0.11936E+01, 0. 122416+01, 0.12726+01 , 0.133496+01, 0.139786+01, 

0.146076+01, 0. 152366*01, 0.158526+01, 0. 164566+01, 0.17036+01 , 

0.17456+01, 0. 178266+01, 0. 181216+01, 0. 182766+01, 0.184066+01 , 

0.1 847 IF + 01 , 0. 18447F + 01, 0.183296+01, 0.181266+01, 0. 178416 + 01, 

C. 1 74836 +01 , 0.169846+01, 0.164096*01, 0.15796E+01, 0.151416+01, 

0. 145986 + 01, 0. 138756 + 01, 0.131526 + 01, 0. 1242 c 5 + 01, 0. 117 066+01, 

0. 1 098 IE + 01 , 0. 10046E + 0 1 , 0.935 86+00, 0.88 286+00, 0.35676+00, 

0. 82766+00, 0.8036+00, 0. 76856+00, 0.73946+00, 0.70096+00, 

0. 6719E+0J, 0.6314E+00, 0.59156+00, 0. 54746+00, 0.50766+00, 

0.4431E + C 0, 0. 39856 + CO , 0.3417E+00, 0.28646+0C , 0.24366+00, 

0. 19956 + 0 J, 0. 13176 + 0O, 0.7056-01, 0.1836-01, -0.3216-01, 

-0.14366+00, -0.15026+00, -0.20386+00, -0.25566+00, -0.29576+00, 

-0.32136+00, -O.3639E+O0, -G.392E+Q0, -O.4334E+O0, -0.47736+00, 

-0.52486 + 00, -0.57386 + 00, -G. 6 198 E+0 v>, -0.6671F + G0, -0.7109E + 00, 

-0.7563F + 00, — 0.7 898E + 00 , -0.83176 + 00, -0.87946 + 00, -0.92776 + 00, 

-0.97936+00, -0.10268E+01, -0.108066*01, -0.113476+01, -0.118236+01, 

-0. 12238E+01, -0.127796+01, -0. 13155E+01, -0. 13722E +01 , -Q. 143936+0 1, 

-0.150126+01, -0.155636+01, -C. 16 26 36 +0 1 , -0.168486+01, -0. 173 68E+0 1 , 

-0.178266+01, -0.181096+01, -0.182636+01, -0. 1 S417E +01, -0. 1 84956+0 1, 

-0. 18466+01, -0.183136+01, -0.180836+01, -0.177886+01, -0.174136+01 , 

-0.169176*01, -0.16446E+01, -0.159186+01, -0. 1 5333E +01, -0.146446+01, 

-G. 137476+01, -0.12983F+01, 
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NKR 

R 

NMAT 
I Ft AG 
W 


KO 

K LOSE 
$ FN D 


$NAK1 

NKP 

P 

$ENC 


= 10, 

* 0.0, 0.1 £*00, 0* 46 + 00 , 0.66+00, Q.8E+00, O.IE+Ol, 0.126*01, 

0.146*01, 0.156*01, 0. 166*01, 0. 171 140736*01 , 

* 157, 

» It 


0.16+01, 

0. 1E+01, 

3. 16+01, 

0. IE +01 , 

0. 1E+01, 

0.16+01, 

O.IE+Ol, 

0. 16+01, 

O.IE+Ol, 

0. 16+01, 

0.1E + 01 , 

0. 16+01, 

0.16+01, 

0.16+01, 

0.16*01, 

0. IE +01, 

0.16+01, 

0.1E+01, 

O.IE+Ol, 

0.16+01, 

0. 1E+01, 

0.16+01, 

0.16+01, 

3. IE +01, 

0.16+01, 

0.16+01, 

0.16 +01, 

0. 1E + 01, 

0. IE + 01 , 

0. 16+01, 

0 . IE + 01 , 

0.16+01, 

0.16+01, 

0.16+01, 

0.16+01, 

J. 16+01, 

0. IE +01, 

0.16+01, 

0.16+01, 

0. 1E+01, 

0. 16+01, 

0.16+01, 

O.'lEfrOi , 

0. lE+Ol, 

0. 16+01 , 

0. IE +01, 

0. IF +01 , 

0.16+01, 

0.16+01, 

0* IE +01 , 

0.16+01, 

0.16+01, 

0. 16+01, 

O.IE+Ol, 

0.16+01, 

3. IE +01, 

0. IE + 01 , 

0. 1E+01, 

0.16*01, 

0. 16+01, 

0. 16 + 01, 

0.16+01, 

0. 16 + 01, 

0. 16+01, 

0. 16 + 01, 

0. 16+01, 

O.IE+Ol, 

0. 16+01, 

0. 16+01, 

0.16+01, 

0.16 +01, 

0. 16+01, 

0. 16 +01 , 

0. 16+01, 

O.IE+Ol, 

0. 16+01, 

0. 16+31, 

0 • IE +0 1 , 

0. IE +01, 

0. IE +01 , 

0. 16 + 01, 

0. IE + 01, 

O.IE+Ol, 

0. IE + 3 1, 

0. 16 +01 , 

0. IE +01 , 

0.16+01, 

0.1E+C1, 

0.16+01, 

0.16+01, 

0. 16+01, 

0.16+01, 

0.16+01, 

0. IE +01 , 

0.16+01, 

0.16+01, 

0 . IE +01, 

3.16+31, 

O.IE+Ol, 

0.16+01, 

0.16 +01, 

0. 16 + 01, 

0. 16 +01, 

0.16+01, 

0.1E+31, 

0. IE *01 , 

0. IE +0 1 , 

J • 16 + 01 , 

0.16+01, 

o. le+ui, 

0.16+01, 

O.IE+Ol, 

0.16+01, 

0.16+01, 

0.16+01, 

0.16+01, 

0. 16+01, 

0.16+01, 

3. If +01, 

0. IE +01, 

0. 1E+01, 

0.16+01, 

0. 16+01 , 

0. 16+01, 

0.16 +01, 

0.16+01, 

0 • IE +01 , 

0.1E+01, 

0 • IE +01, 

0 • 16 +01, 

0. 16 +01, 

0.1E+O1, 

0. 16+01, 

0. 16 + 0.1, 

0. IE + 01, 

0. IE +01, 

0. IE + 01 , 

0. 16+01 , 

0. IE +01, 

0.16+01? 

0. 16+01, 

0.1E+J1 , 

0.16+01, 

0. l'E + 01, 

0.16+01, 

0.16+01, 

0.16+01, 

0. IE +01 , 

0. IE +31, 

0.16+01, 

0. 1E + 01, 

0. 16+01, 

0. 16+01, 

0. 16+01, 

0. IE +01 , 

0. 16+01, 

0. IE +01, 






= 4, 

= 1 , 


The additional input junction values for 


SAy(t) are as follows: 


= to, 

= 0.0, 0.16+00, 0.4E+00, 0.66*00, 0.86+00, 0.16*01, 0.126+01, 

0.146+01, 0.156+01, 0.166*01, 0.17114C136+01, 


The output for this case is as follows: 
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NUMBER OF CURVES = 10 NUMBER OF POINTS = 157 

ENDPOINT S AND JUNCTION POINTS = 


APPENDIX C - Continued 


O 

o 

+ 

Uj 

o 

o 

o 

o 

o 

o 

© 


00 

00 

in 

eg 

eg 

o 

in 


g- 

o 

re 


-4 O 

0 o 

1 + 

UJ UJ 

o re 
o e- 
o o 
o g* 

O -4 
o -4 
O I- 


r- 

O' 

t- 

re 

g- 


+ 

x 

o 

00 

00 

rg 

eg 

eg 


-« o 

rg 

*•4-4 

eg 

rg 

rg 

w4 

-4 

eg eg eg eg re 

eg 

r* 

re 

o o 

o 

o o 

o o 

o 

o 

o 

o 

o o o o o 

o 

00 

re 

i + 

m i 

i i 

i 1 

i 

1 

i 

1 

1 } t 1 1 

i 

• 


UJ uj 

-J UJ 

UJ UJ 

UJ UJ 

UJ 

UJ 

uu 

UJ 

UJ Uj UJ UJ Uj 

UJ 

f-»4 

43 

o o 

< rg 

rg e- 

oo 

o 

rg 

'O 

g 

r- p- in -h cc 

rg 

1 


o o 

O r- 

I s - re 

00 0 s 

*4 

CO 

O' 


cr O' rg e- co 

i"^ 


| 

o-o 

Q 00 

O — » 

in 


CO 

g- 

in • 

• # ^ f 4 in o ^ 

00 



o o 

>4 g* 

g* eg 

g- ^ 

m 

o 

r* 

in 

>f ^oooo 

g 



o © 

oo e- 

o co 

00 v0 

o 

re 

o 

o 

't O' O' vfl 




o.o 

uu o 

rg 

in oo 

O' 

00 

o 

r- 

~4 uvrg <n o 

o 



O g> 

X 00 

O re 

*•4 O' 

o 

CO 

g- 

v0 

r*-. g co ^ 

CO 



• • 

* 

• * 

• • 

* 

• 

* 

• 


t 



vO —< 

m 

—4 —4 

-4 gj 

re 

g- 

<-4 

•H 

in g* m m 

<n 




I I l I t I 


l I 


-4 o 

o 

o 

o 

o 

o 

o o 

o 

o 

-4-4000 

o 

rg 

# 

* 

rg 

* 

* 

o o 

X o 

o 

o 

o 

o 

o O 

o 

o 

o o o o o 

o 

X 

X 

1 + 

♦ 

+ 

♦ 

4- 

+ 

+ + 

4* 

+ 

i i + + + 

.+ 

V 


UJ ULI 

Q UJ 

UJ 

UJ 

LU 

UJ 

UJ UJ 

UJ 

UJ 

UJ UJ uj UJ UJ 

UJ 

e- 

e- 

o o 

UJ rg 

re 

>0 

*4 

00 

43 O' 

O 

^4 

o o n -4 o 

eg 

g 

fH 

o o 

in 

O' 

00 

g 

CO 

g o 

in 

in 

o eg co o g 

in 

o 

re 

o o 

3 eg 

in 

h- 

in 

re 

O' re 

N. 

in • 

• • in g *o O' o' 

rg 

e- 

'i* 

o o 

ex. rg 

n 0 


o 

00 

o -4 

rg 

r- 

4- O' O' O' rg 

eg 


—4 

o © 

3 e: g> 

g 

in 

43 

o 

re re 

g 

O' 

e- r\j rsj o 

g> 

in 

O' 

o o 

O in 

o 

eg 

in 

■— « i 

gj g) 

co 

-4 

co O' m re rg 

in 

*-f 4 

o 

o m 

U rg 

re 

re 

g 

s0 

e- o 

r-4 

re 

O' f- Q -4 eg 

eg 

43 

v0 

* » 

• 

• 

• 

.♦ 

• 

• • 


• 

• 

re 

eg 

g ^ 

*4 

-4 

«-4 

H 

—4 


eg 

eg 

00 O' t-4 -4 -4 

—4 

■* 

• 


. 1 

1 

J 

1 

1 

1 l 



i I I i t 

t 

43 

co 












.-4 

rg 


-4 © 

o 

o 

o 

o 

Q 

o 

o 

o 

O 

O O 

o 

o 

o 

o 

O 

o 

o 

o 

O 

1 + 

+ 

4- 

+ 

+ 

+ 

+ 

4. 

+ 

+ 

UJ UJ 

UJ 

UJ 

UJ 

Uj 

UJ 

UJ 

UJ 

UJ 

UJ 

o o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o o 

X © 

o 

o 

o 

o 

o 

o 

o 

o 

o o 

re 

e- 

o 

in 

K 

o 

re 

eg 

H 

o O 

g- 

gj 

e- 

rg 

O' 

't 

'd* 


eg 

o o 


o 

in 

e- 

CO 

O' 

—4 

't 

in 

o 

re 



n 

>0 

e- 

O' 

o 

r4 

-4 —4 


J 

-4 




H 

rg 

eg 


l I I I I I I I I 


FF>4 

o 

o 

o 

o 




o o 

o 

o 

o 

o 

pH 



( 1 

+ 

+ 

+ 

4- 

o 



UJ Uj 

uu 

Uj 

UJ 

UJ 

1 



O O 

o 

Q 

o 

O 

Uj 



o o 

o 

O 

o 

o 

•o 



1 o o 

o 

o 

o 

o 

cr 



Q O 

e» 

gs 


re 

o 



o 

o 

g* 

in 

g- 

vC 

re 

m 

-H -4 

rg 

-4 

eg 

-4 

rg 

* 

■* 

re 

o 

—4 

rg 

re 

• 

* 


* * 

• 

* 

* 

• 

F mt 

X 

X 

GO O' 

—4 

-4 

—4 

—4 


w 

.w 

1 1 

1 

1 

1 

1 

!! 

o 

o 






Z 

n 

-t 






o 

o 

St 







r- 

go 







rg 

n 






< 

rg 

g- 







ej 

r- 


o 

-g 

rg 

rg 

rg rg 

rg 

-4 

-4 

Q 

3 

3 

3 

O 

o 

o 

o 

3 

o 

O O 

o 

O 

o 

O 

O 

O 

3 

o 

o 

4- 

1 

1 

i 

t i 

i 

1 

1 

4. 

4. 

+ 

4» 

+ 

4- 

UJ 

UJ 

J4 

UJ 

UJ UJ 

UJ 

UJ 

UJ 

ixl 

uU 

UJ 

jj 

UJ 

UJ 

O 

o 

O 

O 

o O 

o 

g* 

W 

O' 

e- 

o 

g- 

-4 

re 

o 


o 

g? 

s0 O' 

g- 

g- 

co 

3 

O 

3 

g* 

r* 

r~ 

o 

e* 

r- 

m 

in g- 

O' 

g* 

O' « 

• • rn 

re 

CD 

g 

gj 

o 

o 

r* 


-4 

m in 

n 


o 

in 

g- 

r- 

re 

s0 

g' 

o 

g* 

g3 

r- 

in in 

n 

3 

00 

—4 

rg 

re 

r~ 

e- 

-4 

o 

—4 

o 

in 

O' g* 

g r 

gj 

O' 

g- 

>n 

o 


O' 

-4 

rg 

g- 

CD 

o 

m r- 

in 

o 

F— 4 

g> 

>0 

g> 

g> 

v0 

r- 

1 r-l 

o -g 


re 

n gj 

co 

f -4 

fH 


F^ 

—4 

<-4 

-4 

—4 


-4 eg re 

g- 

in g> 

e- 

00 

O' 

rg re 

g- in 

gj r- 









in n 

n 

in 

m 

in 









H 

-4 


“4 

t -4 




59 


133.8936163201 X**3>+ -631 . 694313242 ( X**2 ) + 983. 809098345 (XI* -505.926027359 



Y VALUES 


APPENDIX C — Continued 



60 





NUMBER OF CURVES = 10 NUMBER OF POINTS = Ibl 

ENDPOINTS AND JUNCTION POINTS = 

0. 1. 0000000E-01 4 ♦ 0000000 E-01 b.OOOOOOOE-O L 8.0000000E-01 UOOOOOOOE+OO 
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r X y COMPUTED X COMPUTED Y RES X RES Y 


1 0. -1.3 IE + 00 -1.30E+0Q 

2 1.41E-Q2 -1.416+00 -1.226+00 

3 1. 816-02 -I* 46E+00 -1. 196+00 

4 3.666-02 -1.5 7E+00 -1.116+00 

5 5.406-02 -1.696+00 -1.05E+00 

6 6. 756-02 - 1 • 796+00 -1.006+00 

7 8.556-02 -I. 91E+00 -9.436-01 

8 1.06E-01 -2. 04E+00 -8. 836-01 

9 1.206-01 -2.15E+00 -8.36E-01 

10 1. 416—01 -2.296 + 00 -7.82E-01 

11 1.656-01 -2.436+00 -7.286-01 

12 1.836-01 -2.566+00 -6. 856-01 

13 2.01E-01 — 2. 69E+00 -6. 406-01 

14 2.27 6—0 1 -2.846+00 -5.84E-01 

15 2.44E-01 -2.96E+00 -5.356-01 

16 2.626-01 -3.086+00 -4. 826-01 

17 2.70E-01 —3. 16E+00 -4.446-01 

18 2.796-01 -3.246+00 -3.936-01 

19 2.876-01 -3. 306+00 -3.29E-01 

20 2.986-01 -3.346+00 -2.296-01 


-1.266+00 -1.296+00 5.81E-02 

“1.306+00 -1.21E+00 1.02E-01 

-1.336+00 -1.196+00 1.326-01 

-1.466+00 -l.IOE+OO 1.166-01 
-1.626+00 -1.046+00 6.996-02 

-1.766 + 00 - 1 * 0IE + 00 3.096-02 

— 1 . 96E+00 -9. 626-01 4.886-02 

-2.186+00 -9.12E-01 1.406-01 

-2.326+00 -8.756-01 1.68E-01 

-2.516+00 -8.156-01 2 . 19E-G 1 

-2.696+00 -7.386-01 2.546-01 

-2 . 796+0 0-6.806-01 2 . 366-0 1 

-2.896+00 -6.126-01 2.026-01 

— 2* 98E+00 -5.126-01 1.49E-01 

-3. 03 E +00 -4.40E-01 6.93E-02 

—3. 05 £+00 -3.576-01 3.03E-02 

-3.066+00 -3.20E-01 1.066-01 

-3. 066+00 -2.766-01 1.88E-01 

-3.05E+00 -2.40E-01 2.486-01 

-3.04E+00 - 1.836-01 2. 96 6-01 


6.926-03 
1 . 686-02 
1.11E-04 
1.036-02 
9.36E-03 
3.79E-03 
1.856-02 
2.916-02 
3.93E-02 
3. 266-02 
1 . 056-02 
4.87E-03 
2 . 72 E - 02 
7.166-02 
9.586-02 
1 • 256—01 
1.24E-01 
1. 176-01 
8.946-02 
4.59E-02 


147 1.576+00 -2. 77t-01 

148 1.596+00 -4.006-01 

149 1.606+00 -5. 176-01 

150 1.616+00 -6.30E-01 

151 1.63E+00 -7.506-01 

152 1.6464-00 -8. 426-01 

153 1. 656+00 -9.326-01 

154 1.666+00 -1.026+00 

155 1.686+00 -1. 116+00 

156 1.706+00 -1.236+00 

157 1. 716+00 -1.3 1E + 0Q 


-1.83E+00 -2.53E-01 
-1.8 IE +00 -4.036-01 
-1.78E+Q0 -5.41E-Q1 
-1.74E+00 -6.70F-01 
-1.69E+00 -8.146-01 
-1.64E+00 -8.99E-01 
-1.59E+00 -9.79E-01 
-1.53E+00 -1.056+00 
-1.46E+00 -1. L3E+00 
-1.37E+00 -1.22E+00 
-1.306+00 -1. 266+00 


-1.83E+00 2.39E-02 
-1.80E+00 2.75E-03 
-1.76E+00 2.316-02 
-1.73E+00 3.936-02 
-1. 686+00 6.406-02 
-1.646+00 5. 71E-02 
-1.60E+00 4. 76E-02 
-1.556+00 3.43E-02 
-1.49E+00 1.846-02 
-1.38E+00 5. LIE-03 
-1.29E+00 5.816-02 


3.45E-03 
1. 106-02 
1.56E-02 
1.546-02 
1.64E-02 
5. 526-03 
8.286-03 
1.76E-02 
2.336-02 
2.826-03 
6. 926-03 


MAXIMUM RESIDUAL = 3.07040536-01 

A better fit can be realized by changing the joints for S^ x (t) to the following 
values: 


NKR = 20, 

R = 0.0, 0.26+00, 0.256+00, 0.286+00, 0.36+C0, 0.336+00, 0.366+00, 

0.56+0C, 0.76+0C, 0.8756+00, 0.896+00, 0.916+0C, 0.16+01, 

0. 1096*01, 0.1116+01, 0.1136+01, C. 1166+01, 0.1246+01, 

C. 1366+01 , 0.1496+01, 0. 171140736+01 , 


and changing the joints for S^y(t) to the following values: 


NKR * 27, 

R = 0.0, C. 26+00, 0.256*00, 0.286+00, 0.36+00, 0.336+00, 0.386+00, 

0. 446+00, 0 .486+00, 0.526+00, 0.66+00, 0. 76+00, 0.84E+00, 

C. 966 +00 , 0.1046+01, 0.116+01, 0. 1136+01, 0.1166+01, 0.126+01, 

0.1246+01, 0.1326+01, 0.1386+01, C.1446+C1, 0.156+01, 

0.1546+01, 0.1626+01, 0.1696+01, 0.171140736+01, 


The output with these improved spline curves is as follows: 
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32.00 9001117(X**3)+ -151.159273734(X**2)* 229. 9406898791 XI+ -112.542534222 



NUMBER OF CURVES ,=>*4^ 27 ENDPOINTS AND JUNCTION POINTS 



T VALUES 









NUMBER OF CURVES - 27 NUMBER OF POINTS = 157 

ENDPOINTS AND JUNCTION POINTS ~ 

0. 2*00000006-01 2. 5000000E-01 2.80000006-01 3 .OOOOOOQE-Ol 3.30000006-01 

3.80000006-01 4 .40000006-0 1 4.80000006-01 5.20000006-01 6 .00000006-01 7.00000006-01 
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2 47. 3 4 034 77531 X**3) + -12 53 . 767846843 < X**2» *■ 2123. 5 33508 889 ( X) «• -1203.156536949 




X VALUES 
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APPENDIX C - Continued 


T X Y COMPUTED X COMPUTED Y RES X RES Y 

1 0. -1.31E+00 -1.30E+00 -1.336+00 -1.30E+00 1.656-02 3.05E-03 

2 1.41E-02 —1.4 16 + 00 -1.22E+00 -1.426+00 -1.22E+00 1.39E-02 1.57E-03 

3 1.81E-02 — 1. 46E+00 -1.L9E+00 -1.456+00 -1.20E+00 1. 116-02 1.656-02 

4 3.666-0 2 -1.57E+00 -1.116+00 -1.576+00 -1.126+00 4.34E-03 4.32E-03 

5 5.40E-02 — 1. 69E+00 -1.05E+00 -1.69E + 00 -1.05E+00 3. 716-03 5.46E-03 

6 6.75E-02 — 1 . 79E+00 -l.OOE+OO -1.78E+00 -9.99E-01 1.50E-02 3.61E-03 

7 8.556-02 —1.9 IE+00 -9.43E-01 -1.90E+00 -9.39E-01 1.026-02 3.97E-03 

8 1.066-01 -2.04E+00 -8.83E-01 -2.05E+00 -8.78E-01 3.396-03 4.48E-03 

9 1.20E-01 -2.15E+00 -8.36E-01 -2.146+00 -8.40E-01 8.60E-03 4.19E-03 

10 1.41E-01 -2.29E+00 -7.82E-01 -2.29E+00 -7.86E-01 1.95E-03 3.69E-03 

11 1.656—01 —2. 43E+00 -7.28E-01 -2.456+00 -7.27E-01 1.89E-02 6.44E-04 

12 1.83E-G1 -2.56E+00 -6.85E-01 -2. 56E+00 -6.86E-01 5. 786-03 1.67E-03 

13 2. 01 E— 01 -2.69E+00 -6.40E-01 -2.68E+00 -6.42E-01 7.94E-03 2.07E-03 

14 2.27E-01 —2. 84E+00 -5.84E-01 -2.83E+00 -5.786-01 3.10E-03 6.01E-03 

15 2.44E-01 -2.96E+00 -5.35E-01 -2.95E+00 -5.34E-01 8.58E-03 1.13E-03 

16 2.62E-01 -3.086+00 -4.82E-01 -3. 10E+00 -4.846-01 1. 446-02 2.56E-03 

17 2.706-01 -3.16E+00 -4.44E-01 -3. 16E+00 -4.526-01 9.146-04 7. 976-03 

18 2.79E-01 -3.24E+00 -3.93E-01 -3.246+00 -3.94E-01 3.256-03 9.70E-04 

19 2.87E-01 -3.306+00 -3.296-01 -3.30E+00 -3.286-01 5.596-03 1.08E-03 

20 2.986-01 —3 . 34E+00 -2.296-01 -3.346+0 0 -2.066-01 2.466-03 2.28E-02 


147 1.57E+00 -2.776-01 -1.836+00 -2.74E-01 -1.846+00 3.266-03 3.886-03 

148 1.596+00 -4.006-01 -1.816+00 -4.006-01 -1.81E+00 3.16E-04 4.41E-03 

149 1.606+00 —5. 17E-01 -1.78E+00 -5.17E-01 -1.78E+00 7.356-04 2.936-03 

150 1.616+00 -6.306-01 -1.746+00 -6. 296-01 -1.74E+0Q 1.386-03 9.19E-04 

151 1.636+00 -7.50E-01 -1.69E+00 -7.60E-01 -1.68E+00 1.01E-02 8.52E-03 

152 1. 646+00 -8.42E-01 -1.646+00 -8.42E-01 -1.64E+00 8.976-06 4.83E-03 

153 1.656+00 -9.32E-01 -1.596+00 -9.246-01 -1.596+00 7.946-03 4.476-04 

154 1.666+00 -1.02E+00 -1.53E+00 -1.016+00 -1.54E+00 1.34E-02 4.476-03 

155 1.686+00 -1.11E+00 -1.466+00 -1.106+00 -1.476+00 1.036-02 6.32E-03 

156 1.70E+00 -1.236+00 -1.37E+00 -1.24E+00 -1.376+00 1.73E-02 7.176-03 

157 1.71E+00 -1.316+00 -1.306+00 -1.336+00 -1.306+00 1.656-02 3.05E-03 


MAXIMUM RESIDUAL = 5.41674436-02 
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Case 3 

Case 3 is a parametric algorithm, arbitrary curve. The input with junction points 
for S/\ x (t) is as follows: 


tNAMl 

X = -0.1263E>01, - 0. 1 22072E+01 , -0.117843E*-01 , -0. 1 1361 4E+01 , 

-0. 109383E ♦ 01 , -0. 1051 5E + 01 , -0. 1039 15E>3 1 » -0.96676E+00, 
-0.52434E*-00, -0. 88189 E *-00, -0. 83938E+00, -0. 79683E +00, -0. 75423EOO, 
-0.71156E+G0, -0.66883E>00, ~0.62603E>00, -0.583 16E <-00, -0. 54021E*-30» 
-0.497 1 9E *-00 , -0. 45412E*-30» -C.41104E»-00, -0.36795E +00 » -0.32439E*00, 
-0.28189E+03, -0.23896E*-00, -0.19613E+00, -0. 15343E *-00 , -0.11 088E*-0 0, 
-0.6851E-01, -0. 2634E-01 » 0.1561E-01, 0.5731E-01, 0.9373E-01. 

0. 13986E + 00, 0. 18067E+00, 0. 22124E+00, 0. 26 166E *-00, 0.30202E*-90, 
0. 34242E *- 00* 0. 38295E *-00, 0.4237E*-00, 0.46477E*- 00 , 0 . 50625E +00 , 

0. 54324E 4-30, 0. 59082E*-00, 3.6341E>00, 0.67816E*- 00, C.7231E+00, 

0.76901E4-00, 0. 81599E*-00, 0. 86401 E*-00 , 0.90429E *-00, 0. 91447E+00, 

0. 89488E*-00, 3. 87221 E *-00, 0. 35497E*-00, 0. 84167E+00, 0. 32 972E*-00, 

0. 81t 2E + 0 0 , 0. 79989E *-00 , 0.78051E+00, 0.75792E+00 » 0.73271E+00, 

0.70571E*-00, 0. 67775 E*-00, 0.6 4967E *-0 0, 0. 62229E+00, 0.59644E4-00, 
0. 57297E *-00, 0. 55253E*-00, 0. 53474E *-00, 0.5188E*-90». 0.53396F*-09, 

0. 48942E *-00 , 0. 47442E*-00, 3.45816E*-00, 0. 43989E >90» 0.41 87 8 E*C 0, 

0.39439E+00, 0. 37277E *-00, 0.36649E*-Q0, 0. 33424E *-00, 0.41364E+00, 

0.43579E *-00, 0. 43654E*-00, 0.41777E*-00, 0.3847E*-00» 0. 342 52E *-00 , 

C. 29t>46E *• 00, 0.25107E*-00, 0.2036E + 00, 0. 17075E + 00, 0 . 1 3925E *-00 , 

3. 1147EOO, 0.94E-01, 

YY = 0.0, -0.19E-03, -0.36E-03, -0.5E-03, -0.61E-03, -0.65E-03, 

-0.62E-03, -0.51E-03, -0.3E-03, 0.3E-04, 0.49E-03, 0.109E-02, 

0. 185E-02 , 0. 278E-02 , C.389E-02, 0.521E-02, 0. 674E-02, 

0. 85E-02 , 0. 1046E-01, 0.1256E-01, 3.1471E-01, 0.1636E-01, 

0.1893E-01, 0.2085E-01, 0.2254E-01, 0.2393E-01, 0.2496E-01, 

0.2554E-01, 0.2562E-01, 0.251 IE-01, 0.2395E-01, 0.2206E-01, 

0. 1937E-0 1 , 0, 1581F-01 , 0.1135E-01, 0.618E-02, 0.56E-03, 

-C.521E-02 , -0.1 389E-01, -3.1618E-01, -3.2082E-01, -0.2454E-01, 

— 0. 2706E-0 1 , -0. 28 llE-01 » -0.2741E-01, -0.247E-01, -0.197E-01, 
-0.1214E-01, -0.175E-02, 0.1176E-01, 0.2867E-01, 0.5103E-C1, 

0.8355E-01, 3.1 242 7E +00, 0. 16399E*-00 , 0.2014EO0, 0.23937E+00, 
0.28074E + 00, 0.3274E*-00, 0. 37205E+00, 0. 40202 E*- 00 , 0.40676E+00, 
0. 33312E «• CO, 0.35237E *-00, 0. 30573E+00, 0.2546E*-00, 0.205 12E*-00, 
0. 16359E + 00, 0. 13628E *-00, 0. 12845E + 00, 0. 13925E *-30, 0. 16554E*-00» 

3. 2041 8E *-03, 0. 25203E *-00, 0. 3059 5E *-00, 0.3628E+00, 0 .41944E+03 , 
9. 47271E *-00 , 0. 51962E+00, 0.55997E + 00, 0.59635E+00, 0.6313E*00, 

3.66623E *-00, 0. 70223E *-0 0 , 0. 73973E *-00, 0. 77696E +00 , 0. 81168E*-00, 

0.84165E+00, 0. 36463E + 00, 0. 37934E*-00, 0. 391 69E *-00 , 3.90577E*-00, 

0.9277E*-00, 0. 9606 8E *-00 , 0.1E + 01, 

NKR = 9, 

R = 0.0, 0.2E-01, 0.4E-01 , 3.6E-01, 0.8E-01, 0.1E+00, 0. 112E*-00, 

0. 14E *-00 , 0.16F*-00, 0. 1797541E + 00, 

NMAT = 95, 
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I FLAG 
W 


KG 

KLOSE 
$ END 

$NAM1 

NKR 

R 

SEND 


* -40306 863533 16 885 97# 


— 0*589442 1 66 14664 E- 07 1 
0 . 2464 50 8 493 8 5 79- L 0 4-, 

-0*2649628811 0192- 168 
-0.91865312533783+137 
0.13757355373905-275 
0.4910452888251-219, 
0*52424271680201-219 
0.54084143079047-219 
-0.59272771 5756 79 E- 07 
0*11730037971322-275 
0.13228490334101-275 
0.13186622350189-275 
-0.904561 947801 82E-1 2 
—0. 45119336581121 E- 12 
-0.9 04548 0700 140 IE- 12 
-0.45119336581 121E- 12 
-0.90453419222 62E-12, 
-0.45119336 58 1121 E- 12 
-0.90451684499 144E- 12. 

— 0. 451 1 93 3 65 31 121 E- 12 
-0.15521147254684+145 
-0. 38071476425322-125 

— 0.90450 29 67203 64E-1 2 
-0.45119336581 121 E-l 2 
-0. 90447174218 107E- 12 
—0.45119336581 121 E- 12 
-0.904461 3 33 84 02 IE- 12 
-0.4511933 65 31 121E- 12 
-0*90444745605 24E- 12, 
-0.45119336581121 E-l 2 
-0.15412913210712+145 

0. 10996484795022-123 


0.11 5273067 19054-275 
-0. 17365204637693-163 
0.25929515873337-229 
-0.13142742952639-200 
0.471630122 50524-219 
0. 49934464581 933—219 , 
0.51370981213919-219 
U. 54914078778469-219 
0. 1 15273067 19054-2 75 
- 0. 9045792950365 8E— 1 2 
-0.451 1933 65 81 12 IE-12 
-0. 15521444967826+1 45 
-0.3807147642 5322-12 5 
-0.60630411753296+142 
-0.38071476425322-125 
-G. 23683663711337+140 
-0.38071476425322-125, 
-0.9251395637 3211 + 137 
-0.38071476425322-125 
-0.3 97342 39442937+147 
0.11730037971322-275 
0. 1 3229490334 10 1—275 
-0.38071476425322-125 
-0. 6062924981 1337+142 
-0.38071476425322-125 
-0.23683209437134+140 
-0 o 3 807147 642 532 2-125 
-0.92512181864607+137 
-0.38071476425322-125, 
-0.3 9733527297344+147 
0.76435957915019+299 
- 0. 59 1966013 94405E-07 


-0. 23553149769543+140, 
-0.35884981876725+135, 
- 0. 4043 010315618-173, 
0.20445924810045-231 , 
0.48274593183088-219, 
0.50764400231356-219, 
0.53254207379624-219, 
0. 54914078778469-219 , 
-0.397359 13734271 +147 , 
-0. 38071476425322-125 , 
— 0. 59 27 2 3 16 82 8328 E- 07 , 
0. 11730037971322-275, 
0. 13228490334101-275 , 
0*11730037971322-275 , 
0. 13228490334101-275 , 
0. 11730037971322-275 , 
0.13228490334101-275, 
0. 11730037971322-2 75, 
0.13228490334101-275 , 
0.11730037971322-275 , 
— 0* 904485 61 996888E-12 , 
— 0. 45 1 1933658112 1 E-l 2 , 
0.13228490334101-275, 
0. 11730037971322-275 , 
0. 13228490334101-275, 
0. 11730037971322-275, 
0. 13228490334101-275 , 
0.11730037971322-275, 
0.13228490334101-275, 
0. 11730037971322-275 , 
-0.35684771720364+135 , 


4, 

0 , 


The additional input junction values for 


SAy(t) 


are as follows: 


= 9, 

= 0.0, 0.2E-01, 0.4E-01, 0.6E-01, 0.8E-01, G.1E+00, 0.112FO0, 

0.14E*00, 0.16E+00, 0. 1 797541 E>00, 


The output for this case is as follows: 
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NUMBER Of CURVES = 9 NUMBER OF POINTS = 95 

ENDPOINTS AND JUNCTION POINTS = 

0. 2 . 000000QE— 02 4- 0000000E-02 6.00000006-02 8 .00000006-02 1.00000006-01 
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APPENDIX C - Continued 



T 

X 

Y COMPUTED X 

COMPUTED Y 

RES X 

RES Y 

1 

0. 

-1. 266+00 

0. 

-1.266 + 00 

2. 3 IE- 03 

8.47 E— 04 

2.316-03 

2 

1 • 796—03 

-1.22E+Q0 

-1.90E-04 

— 1 . 22 E+00 

-1.25E-04 

2.74E-05 

6.456-05 

3 

3.588-03 

-1. 18E+00 

-3.606-04 

-1.18E+00 

-1. 686-0 3 

4.85 E— 04 

1.32E-03 

4 

5.36E-03 

-1. 14E+00 

-5.00E-04 

-1.14E+00 

— 2.48 E— 03 

7. 22 £-04 

1.986-03 

5 

7.156-03 

-1. G9E+Q0 

-6. 106-04 

-1.09E+00 

-2 . 646-03 

7.476-04 

2.036-03 

6 

8.956-03 

-1.05E+00 

-6.50E-04 

-1.056+00 

-2.27E-03 

6.016-04 

1.626-03 

7 

1.076-02 

-1.01000 

-6.20E-04 

-l. 016+00 

-1. 506-03 

3. 286-04 

8.76E-04 

8 

1.256-02 

-9.676-01 

-5. L0E-04 

-9.676-01 

-4. 35E-04 

1.156-05 

7.496-05 

9 

1.43E-02 

-9.246-01 

-3.00E-04 

-9.25E-01 

7.93E-04 

3.82E-04 

1.096-03 

10 

1.61E-02 

-8.826-01 

3.00E-05 

-8.83E-01 

2.07E-03 

7.41 E— 04 

2.046-03 

11 

l • 79 E— 02 

-8.396-01 

4.90E-04 

-8.406-01 

3.266-03 

1.01E-03 

2.776-03 

12 

1.98E-02 

-7.976-01 

1.09E-03 

-7.98E-01 

4.266-03 

1.17E-03 

3 . 17E-03 

13 

2. 16 E -02 

-7. 546-01 

1.856-03 

-7.556-01 

4. 966-03 

1.176-03 

3. 116-03 

14 

2. 3 4 E- 02 

-7. 126-01 

2.78E-03 

-7. 136-01 

5.396-03 

1 • 00 £-03 

2.61E-03 

15 

2,526-02 

-6.696-01 

3. 896-03 

-6.70E-01 

5.666-03 

7.036-04 

1.776-03 

16 

2.71E-Q2 

-6.26E-0I 

5,216-03 

-6.26E— 01 

5.846-03 

2.97 £—04 

6. 346-04 

17 

2.89E-02 

-5. 836-01 

6.74E-03 

-5. 836-01 

6.056-03 

1.866-04 

6.896-04 

18 

3. 07 £-02 

-5.406-01 

8.506-03 

-5.39E-01 

6.386-03 

7.266-04 

2.126-03 

19 

3.26E-02 

-4.976-01 

1.056-02 

-4.96E-01 

6.926-03 

1.286-03 

3.54E-03 

20 

3. 45 E -02 

• 

—4. 546- 01 

• 

1.266-02 

.# 

-4.52E-01 

• 

7.78E-03 

• 

i, 77E-03 

• 

4.78E-03 

• 

8 5 

1. 59E-01 

• 

4. 376-01 

• 

• 

7.40E-01 

• 

3.926-01 

. ♦ 

7.72E-01 

m 

4.446-02 

• 

3.196-02 

86 

1 .6 IE— 01 

4. 18E-01 

7.T7E-01 

3. 8 IE- 01 

7.98E-01 

3.71E-02 

2.106-02 

8 7 

1.63E-01 

3.856-01 

8. 126-01 

3. 6 IE— 01 

8.226-01 

2.346-02 

9. 926-03 

88 

1.666-01 

3.436-01 

8.426-01 

3. 326-01 

8.39E-01 

1.016-02 

3.026-03 

89 

1. 69E-01 

2. 96E-01 

8.65E-01 

2.97E-01 

8. 516— -01 

4.076-04 

1.346-02 

90 

1. 71E-01 

2. 516-01 

8.806-01 

2.60E-01 

8.64E-01 

9.166-03 

1.576-02 

91 

1. 736-01 

2.09E-01 

8.92E-01 

2. 25 E- 01 

S.8GE-01 

1.616-02 

1. 156-02 

92 

1.75E-01 

1. 7 IE- 01 

9. 066-01 

1.92E-01 

8.996-01 

2.086-02 

6.296-03 

93 

1.76E-01 

1.396-01 

9.28E-01 

1.59E-01 

9.23E-01 

1.976-02 

4.806-03 

94 

l. 786-01 

1. 156-01 

9. 616-01 

1. 18E-01 

9 * 5 8 E— 0 1 

3.66 E-0 3 

2.266-03 

95 

1.80E-01 

9.40E-02 

1 . 00 E+00 

6.656-02 

1,01 E+00 

2.756-02 

1.406-02 

PAX l MUM Rfc $ i DUAL = 

l • 05690661 

:-oi 





A better fit can be obtained by changing the joints for S^ x (t) to the following 
values: 

NKR 
R 


26, 

0.0, 0.4E-01, 0.6E-01, 0.64E-01, G.74E-01, 0. 86E-01, 0.9F-01, 

0.94E-01, 0.98E-01, 0. 102E+00, 0.106E + 00, 0. 109E +00, 0. U2F+00, 

0. 114E + 00, 0. 12E+03, 0.124E+00, 0. 126E+00, 0.129E+00, 

0.1 32E *00 , 0. 138E+00, 0.148E + 00, 0. 152E + 00, 0. 156E+00, 

0. 16E <-00, 0.164E*00, 0. 172F+00, 0. 1 797541 E + 30 , 


and changing the joints for S^y(t) to the following values: 


NKR = 22, 

R = 0.0, 0.75E-01, 0.BE-01, 0.85E-01, 3.9E-01, 0.94E-01, 0.968-01, 

0.1E + 00, 0. 109E+Q0, 0. 109E+00, 0.111E + 00, 0. 1128+00, 3.118F+00, 

0. 1248+00, 0. 126E+00, 0. 127E+00, 0. 128E+00, 0. 132E+00, 

0. 142E+00, 0.15E+00, 0.16E+00, 0. 17E+00, 0. 1 7 97541 E+00, 


The output with these improved spline curves is as follows: 
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NUMBER OF CURVES = 26 ENDPOINTS AND JUNCTION POINTS 



T VALUES 
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